MODULE rad_carma
  
  USE mem_carma
  
  INTEGER,ALLOCATABLE :: indexi(:)
  INTEGER,ALLOCATABLE :: indexj(:)  
  
  CONTAINS

!kmlnew
  SUBROUTINE radcarma(m1,m2,m3,ia,iz,ja,jz,solfac  &
     ,theta_,pi0_,pp_,rv_,RAIN_,LWL_,IWL_,dn0_,rtp_,fthrd_,rtgt_,f13t_,f23t_ &
     ,glat_,glon_,rshort_,rlong_,albedt_,cosz_,rlongup_ &
     ,mynum,fmapt_,aot_,xland_)
!kmlnew
    ! CATT
    use catt_start, only: CATT ! INTENT(IN)

    USE mem_grid,   ONLY:  centlon,	& !INTENT()
  			   dzm, 	  & !INTENT()
  			   dzt, 	  & !INTENT()
  			   idate1,	  & !INTENT()
  			   imonth1,	  & !INTENT()
  			   itime1,	  & !INTENT()
  			   itopo,	  & !INTENT()
  			   iyear1,	  & !INTENT()
  			   ngrid,	  & !INTENT()
  			   nzp, 	  & !INTENT()
  			   plonn,	  & !INTENT()
  			   time 	    !INTENT()
  
    USE mem_scratch, ONLY: nzpmax	    !INTENT()
    USE mem_radiate, ONLY: lonrad	    !INTENT()
  
    USE rconstants,  ONLY: cp,  	  & !INTENT()
  			   cpor,	  & !INTENT()
  			   p00, 	  & !INTENT()
  			   pi180,	  & !INTENT()
  			   stefan	    !INTENT()
  
    
    USE mem_globrad, ONLY: read_rad_data,rad_data_not_read,raddatfn
    USE mem_aerad, ONLY: ngas,nwave,iprocopio
!kmlnew 
    USE mem_leaf          , only: leaf_g

!srf - new aerosol model
    use aer1_list,       nspecies_aer   =>nspecies &
                        ,spc_alloc_aer  =>spc_alloc 
    use mem_aer1 , only: aer1_g ,AEROSOL

   !LFR->para salvar o estado da memoria
    !LFR  USE rams_rad_state, ONLY: WRITE_radiate_stat
    !LFR->  
  
  
  
    IMPLICIT NONE
  
    INTEGER,INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz,mynum
    REAL,INTENT(IN)    :: solfac
  	 
    REAL,INTENT(IN) :: rtgt_(m2,m3)
    REAL,INTENT(IN) :: f13t_(m2,m3)
    REAL,INTENT(IN) :: f23t_(m2,m3)
    REAL,INTENT(IN) :: glat_(m2,m3)
    REAL,INTENT(IN) :: glon_(m2,m3)
    REAL,INTENT(IN) :: cosz_(m2,m3)
    REAL,INTENT(IN) :: albedt_(m2,m3)
    REAL,INTENT(IN) :: fmapt_(m2,m3)
  
!    REAL,INTENT(IN) :: pm_(m1,m2,m3)   ! particulate material (kg[pm]/kg[air])
    REAL,INTENT(IN) :: theta_(m1,m2,m3)
    REAL,INTENT(IN) :: pi0_(m1,m2,m3)
    REAL,INTENT(IN) :: pp_(m1,m2,m3)
    REAL,INTENT(IN) :: rv_(m1,m2,m3)
!kmlnew
    REAL,INTENT(IN) :: RAIN_(m2,m3)
    REAL,INTENT(IN) :: LWL_(m1,m2,m3)
    REAL,INTENT(IN) :: IWL_(m1,m2,m3)
    REAL,INTENT(IN) :: xland_(m2,m3)
!kmlnew    
    REAL,INTENT(IN) :: dn0_(m1,m2,m3)
    REAL,INTENT(IN) :: rtp_(m1,m2,m3)
  
    REAL,INTENT(INOUT) :: fthrd_(m1,m2,m3)  
    REAL,INTENT(INOUT) :: rshort_(m2,m3)
    REAL,INTENT(INOUT) :: rlong_(m2,m3)
    REAL,INTENT(INOUT) :: rlongup_(m2,m3)
    REAL,INTENT(OUT) :: aot_(m2,m3,nwave)
    
    
    REAL :: rtgt((iz-ia+1)*(jz-ja+1))
    REAL :: f13t((iz-ia+1)*(jz-ja+1))
    REAL :: f23t((iz-ia+1)*(jz-ja+1))
    REAL :: glat((iz-ia+1)*(jz-ja+1))
    REAL :: glon((iz-ia+1)*(jz-ja+1))
    REAL :: cosz((iz-ia+1)*(jz-ja+1))
    REAL :: albedt((iz-ia+1)*(jz-ja+1))
    REAL :: fmapt((iz-ia+1)*(jz-ja+1))
  
    REAL :: pm((iz-ia+1)*(jz-ja+1),m1)! particulate material (kg[pm]/kg[air])
    REAL :: theta((iz-ia+1)*(jz-ja+1),m1)
    REAL :: pi0((iz-ia+1)*(jz-ja+1),m1)
    REAL :: pp((iz-ia+1)*(jz-ja+1),m1)
    REAL :: rv((iz-ia+1)*(jz-ja+1),m1)
!kmlnew
    REAL :: RAIN((iz-ia+1)*(jz-ja+1))
    REAL :: LWL((iz-ia+1)*(jz-ja+1),m1)
    REAL :: IWL((iz-ia+1)*(jz-ja+1),m1)
    REAL :: xland((iz-ia+1)*(jz-ja+1))
!kmlnew    
    REAL :: dn0((iz-ia+1)*(jz-ja+1),m1)
    REAL :: rtp((iz-ia+1)*(jz-ja+1),m1)
  
    REAL :: fthrd((iz-ia+1)*(jz-ja+1),m1)   
    REAL :: rshort((iz-ia+1)*(jz-ja+1))     
    REAL :: rlong((iz-ia+1)*(jz-ja+1))      
    REAL :: rlongup((iz-ia+1)*(jz-ja+1))    
    REAL :: aotl((iz-ia+1)*(jz-ja+1),nwave)
      
    REAL :: prd((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: temprd((iz-ia+1)*(jz-ja+1),nzpmax+1)
    REAL :: dn0r((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: dztr((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: pmr((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: rvr((iz-ia+1)*(jz-ja+1),nzpmax)
!kmlnew
    REAL :: RAINr((iz-ia+1)*(jz-ja+1))
    REAL :: LWLr((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: IWLr((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: xlandr((iz-ia+1)*(jz-ja+1))
!kmlnew
    REAL :: fthrl((iz-ia+1)*(jz-ja+1),nzpmax)
    REAL :: fthrs((iz-ia+1)*(jz-ja+1),nzpmax)
    
    REAL,PARAMETER :: fcui=1.e-9     !de billion g [gas/part] /kg [ar] para kg/kg
    
    
    REAL :: pird,dzsdx,dzsdy,dlon,a1,a2,dayhr,gglon,dztri
    REAL :: dayhrr,hrangl,sinz,sazmut,slazim,slangl,cosi
    INTEGER :: igas,kk,ik,iend,ij,i,j,k,nzz,ispc,imode

    INTEGER :: ncall = 0 

    iend=(iz-ia+1)*(jz-ja+1) !Size of vector
    
    CALL AllocIndex(ia,ja,iz,jz,1) !1 to alloc auxiliar
    
    
    CALL C_2d_1d(rtgt_(ia:iz,ja:jz),rtgt,ia,iz,ja,jz,iend)
    CALL C_2d_1d(f13t_(ia:iz,ja:jz),f13t,ia,iz,ja,jz,iend)
    CALL C_2d_1d(f23t_(ia:iz,ja:jz),f23t,ia,iz,ja,jz,iend)
    CALL C_2d_1d(glat_(ia:iz,ja:jz),glat,ia,iz,ja,jz,iend)
    CALL C_2d_1d(glon_(ia:iz,ja:jz),glon,ia,iz,ja,jz,iend)
    CALL C_2d_1d(cosz_(ia:iz,ja:jz),cosz,ia,iz,ja,jz,iend)
    CALL C_2d_1d(albedt_(ia:iz,ja:jz),albedt,ia,iz,ja,jz,iend)
    CALL C_2d_1d(fmapt_(ia:iz,ja:jz),fmapt,ia,iz,ja,jz,iend)

    if (CATT==1 .and. AEROSOL > 0) then
    !- PM2.5 from biomass burning
       ispc =bburn
       imode=accum
       if(spc_alloc_aer(transport,imode,ispc) == ON) then 
!orig      CALL C_3d_2d(pm_(:,ia:iz,ja:jz),pm,m1,ia,iz,ja,jz,iend)
           CALL C_3d_2d(aer1_g(imode,ispc,ngrid)%sc_p(:,ia:iz,ja:jz),&
	                pm,m1,ia,iz,ja,jz,iend)
       endif
    endif

    CALL C_3d_2d(theta_(:,ia:iz,ja:jz),theta,m1,ia,iz,ja,jz,iend)
    CALL C_3d_2d(pi0_(:,ia:iz,ja:jz),pi0,m1,ia,iz,ja,jz,iend)
    CALL C_3d_2d(pp_(:,ia:iz,ja:jz),pp,m1,ia,iz,ja,jz,iend)
    CALL C_3d_2d(rv_(:,ia:iz,ja:jz),rv,m1,ia,iz,ja,jz,iend)
!kmlnew
    CALL C_3d_2d(LWL_(:,ia:iz,ja:jz),LWL,m1,ia,iz,ja,jz,iend)
    CALL C_3d_2d(IWL_(:,ia:iz,ja:jz),IWL,m1,ia,iz,ja,jz,iend)
    CALL C_2d_1d( RAIN_(ia:iz,ja:jz),RAIN ,ia,iz,ja,jz,iend)
    CALL C_2d_1d(xland_(ia:iz,ja:jz),xland,ia,iz,ja,jz,iend)
!kmlnew

    CALL C_3d_2d(dn0_(:,ia:iz,ja:jz),dn0,m1,ia,iz,ja,jz,iend)
    CALL C_3d_2d(rtp_(:,ia:iz,ja:jz),rtp,m1,ia,iz,ja,jz,iend)
    
    CALL C_2d_1d(rshort_(ia:iz,ja:jz),rshort,ia,iz,ja,jz,iend)       
    CALL C_2d_1d(rlong_(ia:iz,ja:jz),rlong,ia,iz,ja,jz,iend)         
    CALL C_2d_1d(rlongup_(ia:iz,ja:jz),rlongup,ia,iz,ja,jz,iend)     
    CALL C_3d_2d(fthrd_(:,ia:iz,ja:jz),fthrd,m1,ia,iz,ja,jz,iend)    
    CALL Ci_3d_2d(aot_(ia:iz,ja:jz,:),aotl,nwave,ia,iz,ja,jz,iend)
  
      
    ! liga radiacao de onda longa
    ir_aerad = 1
    
!srf ---- otimizacao (comentar o "call end_carma" abaixo)
    !IF (nCALL == 0) THEN
    !  nCALL = 1
    !  CALL init_carma(ia,iz,ja,jz,m1,m2,m3)
    !ENDIF
!    CALL init_carma(ia,iz,ja,jz,m1,m2,m3)
  
    IF(rad_data_not_read) THEN
         CALL init_carma(ia,iz,ja,jz,m1,m2,m3)
  	 CALL read_rad_data
  	 CALL setupbins
    END IF
  
    isl_aerad = 0
    nzz = m1 - 1
    DO k = 1,m1
      DO ij=1,iend
  	pird = (pp(ij,k) + pi0(ij,k)) / cp
  	temprd(ij,k) = theta(ij,k) * pird ! air temperature (K)
  	rvr(ij,k) = max(0.,rv(ij,k))
        
	LWL(ij,k) = max(0.,LWL(ij,k))
	IWL(ij,k) = max(0.,IWL(ij,k))
	
  	! Convert the next 7 variables to cgs for now.
  	prd(ij,k) = pird ** cpor * p00 * 10. ! pressure
  	dn0r(ij,k) = dn0(ij,k) * 1.e-3        ! air density
  	dztr(ij,k) = dzt(k) / rtgt(ij) * 1.e-2

	
      END DO
    END DO


    if (CATT==1 .and. AEROSOL >0) then
       
      DO k = 1,m1
        DO ij=1,iend
            pmr(ij,k) = pm(ij,k)*fcui
	ENDDO
      ENDDO
    else
      DO k = 1,m1
        DO ij=1,iend
           pmr(ij,k) = 0.
	ENDDO
      ENDDO
    endif


    
    DO ij=1,iend
      temprd(ij,1) = (rlongup(ij) / stefan) ** 0.25     
      temprd(ij,nzp+1) = temprd(ij,nzp)
      !  Initialize atmospheric structure. 
      p_surf(ij) = 0.5*(prd(ij,1) + prd(ij,2))    
      p_top(ij)  = prd(ij,m1)		  
      t_surf(ij) = temprd(ij,1)   
    END DO
    DO k=1,m1-1
      DO ij=1,iend
  	! K level in CARMA grid corresponds to K+1 level in BRAMS grid
  	! Transfer values from RAMS grid to CARMA grid
  	p(ij,k)    =	prd(ij,k+1)
  	t(ij,k)    = temprd(ij,k+1)
  	rhoa(ij,k) =   dn0r(ij,k+1)
!kmlnew
	LWLr(ij,k) = LWL(ij,k+1)*dn0r(ij,k+1) * 1.e+3  ![kg/m3]
	IWLr(ij,k) = IWL(ij,k+1)*dn0r(ij,k+1) * 1.e+3  ![kg/m3]
!kmlnew	
      END DO
    END DO
    
 
    DO ik = 1,NZZ
      !  Reverse the vertical index when in cartesian coordinates
       kk = NZZ + 1 - ik
      DO ij=1,iend
  	 t_aerad(ij,kk) = t(ij,ik)
  	 p_aerad(ij,kk) = p(ij,ik)
!kmlnew
	LWL_aerad(ij,kk) = LWLr(ij,ik)
	IWL_aerad(ij,kk) = IWLr(ij,ik)
	
!        if(LWL_aerad(ij,ik).ne.0.) print*,'LWL [g/cm3]=', LWL_aerad(ij,ik), ik
!	if(IWL_aerad(ij,ik).ne.0.) print*,'IWL [g/cm3]=', IWL_aerad(ij,ik), ik
!kmlnew	
	 
      END DO
    END DO

!kmlnew 

    DO ik = 1,NZZ
      DO ij=1,iend
        dztri=1./(dztr(ij,ik) * 1.e+2)
	LWP_aerad(ij,ik) = LWL_aerad(ij,ik) * dztri   ![kg/m2]
	IWP_aerad(ij,ik) = IWL_aerad(ij,ik) * dztri   ![kg/m2]
!	if(IWL_aerad(ij,ik).gt.0.) print*,'ik,IWL=',ik,IWL_aerad(ij,ik),IWP_aerad(ij,ik)
      END DO
    END DO
!kmlnew 
      DO ij=1,iend
       xland_aerad(ij)=xland(ij)
       RAIN_aerad(ij)=RAIN(ij)
      END DO
   
    
    DO ij=1,iend
      tabove_aerad(ij)  = t(ij,nzz)
    END DO
    
    !  Initialize gas concentrations.
    DO igas = 1,ngas
       DO k = 1,m1-1
  	  ! K level in CARMA grid corresponds to K+1 level in BRAMS grid
  	  DO ij=1,iend
  	    ! water vapor concentration
  	    IF( igas .eq. 1 ) gc(ij,k,igas) = rvr(ij,k+1) * dn0r(ij,k+1) 
  	  END DO
       END DO
    END DO
  
    DO ij=1,iend
      ! The shortwave parameterizations are only valid if the cosine
      !    of the zenith angle is greater than .03 .
      IF (cosz(ij) .gt. .03) isl_aerad(ij) = 1
    END DO
      
    CALL initaer(m1,pmr,dn0r,ia,iz,ja,jz,nzpmax)
    
    !  Initialize radiation
    CALL initrad(imonth1,idate1,iyear1,itime1,time,m1,ia,ja,iz,jz)
  
    CALL prerad(m1,dztr,fmapt,ia,iz,ja,jz,nzpmax,m2,m3)
    
    CALL radtran(albedt,cosz,m1,m2,m3,ia,iz,ja,jz,aotl(:,11))  !kml2
  
    CALL radtran_to_rams(nzp,m2,m3,fthrl,rlong,fthrs,rshort,aotl,ia,iz,ja,jz,mynum)   
  
    ! 
    ! Modify the DOwnward surface shortwave flux by considering
    !	 the slope of the topography.
  
    DO ij=1,iend
      IF (itopo .eq. 1) THEN
  	dzsdx = f13t(ij) * rtgt(ij)
  	dzsdy = f23t(ij) * rtgt(ij)
  
  	! The y- and x-directions must be true north and east for
  	! this correction. the following rotates the model y/x
  	! to the true north/east.   
  
  	! The following rotation seems to be incorrect,so CALL this instead:
  	! SUBROUTINE uvtoueve(u,v,ue,ve,qlat,qlon,platn(ngrid),plonn(ngrid))
  
  	dlon = (plonn(ngrid) - glon(ij)) * pi180
  	a1 = dzsdx*cos(dlon) + dzsdy * sin(dlon)
  	a2 = -dzsdx*sin(dlon) + dzsdy * cos(dlon)
  	dzsdx = a1
  	dzsdy = a2
  
  	dayhr = time / 3600. + float(itime1/100)  &
  	    + float(mod(itime1,100)) / 60.
  	gglon = glon(ij)
  	IF (lonrad .eq. 0) gglon = centlon(1)
  	dayhrr = mod(dayhr+gglon/15.+24.,24.)
  	hrangl = 15. * (dayhrr - 12.) * pi180
        !srf - evitando SQRT (<0)
        !sinz = sqrt(1. - cosz(ij) ** 2)
  	sinz = sqrt(max(0., (1. - cosz(ij) ** 2)))
  
  	! ALF - Evitando divisao por zero
  	sinz = max(0.000001, sinz)
  	! ALF
  
  	sazmut = asin(max(-1.,min(1.,cdec*sin(hrangl)/sinz)))
  	IF (abs(dzsdx) .lt. 1e-20) dzsdx = 1.e-20
  	IF (abs(dzsdy) .lt. 1e-20) dzsdy = 1.e-20
  	slazim = 1.571 - atan2(dzsdy,dzsdx)
  	slangl = atan(sqrt(dzsdx*dzsdx+dzsdy*dzsdy))
  	cosi = cos(slangl) * cosz(ij) + sin(slangl) * sinz  &
  	     * cos(sazmut-slazim)
  	rshort(ij) = rshort(ij) * cosi / cosz(ij)        
     END IF
    END DO
    
    !print*,'------ radiative heating rates ---- ---'
    DO k = 2,m1-1
      DO ij=1,iend
  	 fthrd(ij,k) = fthrl(ij,k) + fthrs(ij,k)        
!	 if(fthrd(ij,k)*86400. < -10.) &
!	 print*,'IJ K',ij,k,fthrl(ij,k)*86400. , fthrs(ij,k)*86400.
      END DO
    END DO
  
    ! Convert the downward flux at the ground to SI.
  
    !	     rshort(i,j) = rshort(i,j) * 1.e-3 / (1. - albedt(i,j))
    !	     rlong(i,j) = rlong(i,j) * 1.e-3
    DO ij=1,iend
      rshort(ij) = rshort(ij) / (1. - albedt(ij)) 
      rlong(ij) = rlong(ij)                       
      fthrd(ij,1) = fthrd(ij,2)                  
      !print*,'IJ SW LW=',ij,rshort(ij),rlong(ij)
      !call flush(6)
      !if (rlong(ij).lt.10.) print*,'Antes da conversao!!!',ij,rlong(ij)

    END DO
  
    
!srf ---- comentar a linha abaixo    
    !!! CALL end_carma()
  
    CALL C_1d_2d(rshort,rshort_(ia:iz,ja:jz),ia,iz,ja,jz,iend)    
    CALL C_1d_2d(rlong,rlong_(ia:iz,ja:jz),ia,iz,ja,jz,iend)      
    CALL C_1d_2d(rlongup,rlongup_(ia:iz,ja:jz),ia,iz,ja,jz,iend)  
    CALL C_2d_3d(fthrd,fthrd_(:,ia:iz,ja:jz),m1,ia,iz,ja,jz,iend) 
    CALL Ci_2d_3d(aotl,aot_(ia:iz,ja:jz,:),nwave,ia,iz,ja,jz,iend)
  
    CALL AllocIndex(ia,ja,iz,jz,0) !0 to dealloc auxiliar
    
   
  
  END SUBROUTINE radcarma
  
  SUBROUTINE setupbins
    !  This routine evaluates the derived mapping arrays and sets up
    !  the particle size bins.
  
    USE mem_aerad, ONLY: lunoprt, nbin
    USE mem_globaer, ONLY: ngroup,nelem,itype,i_involatile, &
  			   i_volatile,ienconc,igelem,ncore, &
  			   nelemg,i_coremass,i_volcore, &
  			   i_core2mom,ixyz,nxyz,rhop3,rhoelem,rhopcore3, &
  			   rhocore,pi,rmassmin,rmin,rmrat,one,rmass, &
  			   rmasscore,pcore,rmassup,rmasscoreup,dm,vol, &
  			   r,rcore,rup,rcoreup,dr,rlow,diffmass
    
    IMPLICIT NONE
    
    !Local
  
    INTEGER :: igrp
    INTEGER :: ielem
    INTEGER :: j
    INTEGER :: ie
    INTEGER :: ig
    INTEGER :: ibin
    REAL    :: cpi
    REAL    :: vrfact
  
    !  Determine which elements are particle number concentrations
    !  <ienconc(igroup)> is the element corresponding to particle number 
    !  concentration in group <igroup>
    !
    igrp = 0
    DO ielem = 1, NELEM
       IF( itype(ielem) .eq. I_INVOLATILE .or. &
  	    itype(ielem) .eq. I_VOLATILE )THEN
  
  	  igrp = igrp + 1
  	  ienconc(igrp) = ielem
       END IF
    END DO
    !
    !  Determine which group each element belongs to
    !  i.e., <igelem(ielem)> is the group to which element <ielem> belongs
    !
    igrp = 0
    DO ielem = 1, NELEM
       IF( itype(ielem) .eq. I_INVOLATILE .or.       &
  	    itype(ielem) .eq. I_VOLATILE )THEN
  	  igrp = igrp + 1
       END IF
       igelem(ielem) = igrp
    END DO
    !
    !  Particle mass densities (NXYZ*NBIN for each group) -- the user might want
    !  to modIFy this (this code segment DOes not appear in setupaer SUBROUTINE
    !  because <igelem> is not defined until this SUBROUTINE).
    !
    DO ie = 1,NELEM
       ig = igelem(ie)
       DO ibin = 1,NBIN
  	  DO ixyz = 1,NXYZ
  	     rhop3(ixyz,ibin,ig) = rhoelem(ie)
  	     rhopcore3(ixyz,ibin,ig) = rhocore(ie)
  	  END DO
       END DO
    END DO
    !
    !
    !  Set up the particle bins.
    !  For each particle group, the mass of a particle in
    !  bin j is <rmrat> times that in bin j-1
    !
    !	 rmass(NBIN,NGROUP)	=  bin center mass [g]
    !	 r(NBIN,NGROUP) 	=  bin mean (volume-weighted) radius [cm]
    !	 vol(NBIN,NGROUP)	=  bin center volume [cm^3]
    !	 dr(NBIN,NGROUP)	=  bin width in radius space [cm]
    !	 dv(NBIN,NGROUP)	=  bin width in volume space [cm^3]
    !	 dm(NBIN,NGROUP)	=  bin width in mass space [g]
    !
    cpi = 4./3.*PI
  
    DO igrp = 1, NGROUP
  
       rmassmin(igrp) = cpi*rhop3(1,1,igrp)*rmin(igrp)**3
  
       vrfact = ( (3./2./PI/(rmrat(igrp)+1.))**(ONE/3.) )*    &
  	    ( rmrat(igrp)**(ONE/3.) - 1. )
  
       DO j = 1, NBIN
  	  !PRINT *,'LFRDBG->igrp,j,rmassmin(igrp),rmrat(igrp):', &
  	  !	    igrp,j,rmassmin(igrp),rmrat(igrp)
  	  !CALL FLUSH(6)
  	  rmass(j,igrp)   = rmassmin(igrp) * rmrat(igrp)**(j-1)
  	  rmasscore(j,igrp) = pcore/100. * rmass(j,igrp)
  
  	  rmassup(j,igrp) = 2.*rmrat(igrp)/(rmrat(igrp)+1.)*rmass(j,igrp)
  	  rmasscoreup(j,igrp) = pcore/100. * rmassup(j,igrp)
  
  	  dm(j,igrp)	  = 2.*(rmrat(igrp)-1.)/(rmrat(igrp)+1.)*rmass(j,igrp)
  	  vol(j,igrp) = rmass(j,igrp) / rhop3(1,1,igrp)
  
  	  r(j,igrp)	  = ( rmass(j,igrp)/rhop3(1,1,igrp)/cpi )**(ONE/3.)
  	  !PRINT *,'LFRDBG->cpi,ONE,rmasscore(j,igrp),rhopcore3(1,1,igrp):', &
  	  !	    cpi,one,rmasscore(j,igrp),rhopcore3(1,1,igrp)
  	  !CALL FLUSH(6)
  	  rcore(j,igrp)   = ( rmasscore(j,igrp)/rhopcore3(1,1,igrp)/cpi )**(ONE/3.)
  
  	  rup(j,igrp)	  = ( rmassup(j,igrp)/rhop3(1,1,igrp)/cpi )**(ONE/3.)
  	  rcoreup(j,igrp) = ( rmasscoreup(j,igrp)/rhopcore3(1,1,igrp)/cpi )**(ONE/3.)
  	  !PRINT *,'LFRDBG->vrfact,rmass(j,igrp),rhop3(1,1,igrp):', &
  	  !		    vrfact,rmass(j,igrp),rhop3(1,1,igrp)
  
  	  dr(j,igrp)  = vrfact*(rmass(j,igrp)/rhop3(1,1,igrp))**(ONE/3.)
  	  rlow(j,igrp) = rup(j,igrp) - dr(j,igrp)
       END DO
    END DO
    !PRINT *,'LFRDBG->End of setupbins';CALL FLUSH(6)
  
  END SUBROUTINE setupbins
  
  SUBROUTINE initaer(m1,pmr,dn0r,ia,iz,ja,jz,nzpmax)

    use mem_aerad, only: &
         nbin

    USE mem_globaer, ONLY: nelem      , &
  			   igelem     , &
  			   ienconc    , &
  			   small_pc   , &
  			   itype      , &
  			   i_coremass , &
  			   rmass      , &
  			   fix_coref  , &
  			   i_core2mom , &
  			   rhop3      , &
  			   pi	      , &
  			   dr	      , &
  			   r
    
    IMPLICIT NONE
  
    INTEGER,INTENT(IN)  	      :: m1,ia,iz,ja,jz,nzpmax
    REAL   ,INTENT(IN), DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax) :: pmr
    REAL   ,INTENT(IN), DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax) :: dn0r
  
    !Local
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),m1) :: totm
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),m1) :: r0
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),m1) :: rsig
    INTEGER :: ie,ix,iy
    INTEGER :: ielem
    INTEGER :: ig
    INTEGER :: ip,iend,ij
    INTEGER :: j
    INTEGER :: k
    INTEGER :: kr
    INTEGER :: nzz
    REAL    :: arg1
    REAL    :: arg2
    REAL    :: sum((iz-ia+1)*(jz-ja+1))
    REAL    :: totn((iz-ia+1)*(jz-ja+1))
  
    iend=(iz-ia+1)*(jz-ja+1)
      
    nzz = m1 - 1
    !
    !transfere valores da grade do rams para carma
    DO k = 1,nzz
      kr = K + 1     ! nivel K da grade do carma orresponde ao nivel K + 1 DO RAMS
      DO ij=1,iend
  	!    totm = total mass particle concentration (g/cm3) 
  	totm(ij,k)  = pmr(ij,kr) *  dn0r(ij,kr) 
      END DO
    END DO
    !  Initialize particle number densities 
    !  Core mass is assumed to be 100% of particle mass
    DO ielem = 1,nelem
       ig = igelem(ielem)
       ip = ienconc(ig)
       DO j = 1,nbin
  	  DO k = 1,nzz
  	     DO ij = 1,iend
  	       IF( ielem .eq. ip )THEN
  		  !  Particle number concentration [#/cm^3]
  		  pc(ij,k,j,ielem) = SMALL_PC
  	       ELSE IF( itype(ielem) .eq. I_COREMASS )THEN
  		  !  Core mass concentration [g/cm^3]
  		  pc(ij,k,j,ielem) = pc(ij,k,j,ip)*rmass(j,ig) * &
  					FIX_COREF
  	       ELSE IF( itype(ielem) .eq. I_CORE2MOM )THEN
  		  !  Second moment of core mass distribution [ (g/cm^3)^2 ]
  		  pc(ij,k,j,ielem) = pc(ij,k,j,ip) *	      &
  		       (rmass(j,ig)*FIX_COREF)**2
  	       END IF
  	    END DO
  	  END DO
       END DO
    END DO
    !
    !  Initial particle distribution: log-normal size distribution 
    !  for first particle group (which has only one particle element)
    !  in a single column
    ig = 1
    ie = ienconc(ig)
  
    DO k = 1,nzz
       !
       !  Log-normal parameters:
       !  
       !    r0   = number mode radius
       !    rsig = geometric standard deviation
       !    totm = total mass particle concentration (g/cm3) (proveniente DO rams) 
       !
       DO ij=1,iend
  	 r0(ij,k)  = 1.95e-5
  	 rsig(ij,k) = 1.62
  	 totn(ij) = (6. * totm(ij,k)/(rhop3(1,1,1)*PI*r0(ij,k)**3))* &
  	    exp((-9./2)*log(rsig(ij,k))**2)
       END DO
       !  Adjust prefactor to yield particle number concentration <ntot>
       !
     sum = 0.
     DO j = 1,nbin
        DO ij=1,iend
  	  arg1 = dr(j,ig) / ( sqrt(2.*PI) * r(j,ig) * log(rsig(ij,k)) ) 
  	  arg2 = -log( r(j,ig) / r0(ij,k) )**2 / &
    		       ( 2.*log(rsig(ij,k))**2 )

  	  sum(ij)  = sum(ij) + arg1 * exp( arg2 )
        END DO
     END DO
     DO ij=1,iend
       totn(ij) = totn(ij) / sum(ij)
     END DO
     DO j = 1,nbin
        DO ij=1,iend
  	  arg1 = totn(ij) * dr(j,ig) / ( sqrt(2.*PI) * r(j,ig) * &
  		       log(rsig(ij,k)) ) 

  	  arg2 = -log( r(j,ig) / r0(ij,k) )**2 / &
  		       ( 2.*log(rsig(ij,k))**2 )
  	  pc(ij,k,j,ie) = max( arg1 * exp( arg2 ), REAL(SMALL_PC) )
        END DO
     END DO
    END DO
    
  END SUBROUTINE initaer
  
  SUBROUTINE initrad(imonth1,idate1,iyear1,itime1,time_rams,m1,ia,ja,iz,jz)
  
    USE mem_aerad, ONLY: is_grp_ice_aerad,r_aerad, &
  			 rup_aerad,rcore_aerad,rcoreup_aerad, &
  			 ptop_aerad,pbot_aerad,u0_aerad, &
  			 sfc_alb_aerad,emisir_aerad,tsfc_aerad, &
  			 tabove_aerad,wave_aerad,iprocopio&
                         ,nx,ny,nbin,nwave,nsol
  
    USE mem_globaer, ONLY: time,do_solar,do_ir,isolar_zen,i_diurnal, &
  			   rad_start,scday,pi,ix,iy,rlat,u0,ngroup, &
  			   is_grp_ice,ienconc,r,rup,rcore,rcoreup, &
  			   t_surf,wave,z_sin, &
  			   z_cos
    USE mem_globrad, ONLY: imie
    
    USE mem_carma, ONLY: declin
    
    IMPLICIT NONE
    
    INTEGER,INTENT(IN) :: m1,ia,ja,iz,jz
    INTEGER,INTENT(IN) :: imonth1
    INTEGER,INTENT(IN) :: idate1
    INTEGER,INTENT(IN) :: iyear1
    INTEGER,INTENT(IN) :: itime1
    REAL   ,INTENT(IN) :: time_rams
  
    !Local  
    INTEGER :: iday
    INTEGER :: iwave
    INTEGER :: julday
    REAL    :: saz
    REAL    :: wavetemp
   
    !
    !  Define flag to control the calculation of the solar zenith angle:
    !	 <isolar_zen> = I_FIXED: USE fixed value <u0fixed>
    !		      = I_DIURNAL: calculation based on time, day, lat, and lon
    !
    isolar_zen = I_DIURNAL
  
    IF( isolar_zen .eq. I_DIURNAL )THEN
       !
       !
       !  Define values needed for calculation of solar zenith angle:
       !    <iday> is day of year
       !    <rad_start> is solar time corresponding to <time> = 0, in seconds
       !     = 0 means <time> = 0 corresponds to midnight,
       !     = 6 * 3600 means <time> = 0 corresponds to 6 AM
       !    Note: all times are local standard time.
       !
  
       iday =  julday(imonth1,idate1,iyear1)
       iday =  iday + nint(time_rams/86400.)
       rad_start = (float(itime1/100) + float(mod(itime1,100)) / 60.)*scday
       !
       !
       !  Precalculate terms in solar zenith angle computation:
       !    (adapted from original Toon model)
       !    <saz> is solar azimuth angle [rad]
       !    <declin> is solar declination [rad]
       !    <z_sin> is sin term in precalculation
       !    <z_cos> is cos term in precalculation
       !
  
       saz = 2. * PI / 365. * iday 
  
       declin = 0.006918 - 0.399912*cos(saz)	+0.070257*sin(saz)    &
  	    - 0.006758*cos(2.*saz) +0.000907*sin(2.*saz) &
  	    - 0.002697*cos(3.*saz) +0.001480*sin(3.*saz)
  
       !DO ix = 1,NX
       !   DO iy = 1,NY
       !      rlat(ix,iy) = glat
       !      z_sin(ix,iy) = sin(declin) * sin( rlat(ix,iy) * PI/180. )
       !      z_cos(ix,iy) = cos(declin) * cos( rlat(ix,iy) * PI/180. )
       !   END DO
       !END DO
  
    END IF
    !
  
    !
    !  Initialize the radiative transfer model
    CALL setuprad(m1,ia,ja,iz,jz)
    
     IF(imie == 0) THEN
      CALL calcproperties
     END IF

    !  Get radiative wavelengths
    !
    DO iwave = 1,NWAVE
       !
       !
       !  Solar wavelengths in radiative transfer model are bin centers,
       !  infrared are bin edges
       !
       IF( iwave .le. NSOL )THEN
  	  wave(iwave) = wave_aerad(iwave)
       ELSE
  	  wave(iwave) = 0.5*( wave_aerad(iwave) + wave_aerad(iwave+1) )
       END IF
  
    END DO
    !
    !  Switch bins 11 and 12
    !KLF Corrigidos os comp. de onda (11 e 12) para (17 e 18)!!!
    !	   wavetemp = wave(11)
    !	   wave(11) = wave(12)
    !	   wave(12) = wavetemp
  
    wavetemp = wave(17)
    wave(17) = wave(18)
    wave(18) = wavetemp
    !
    !  Return to CALLer with radiation model initialized
    !
  END SUBROUTINE initrad
  
  SUBROUTINE setuprad(m1,ia,ja,iz,jz)
    !	  *********************************************************
    !	  *  Purpose		:  Defines all constants, and	  *
    !	  *			   calculates pressure averaged   *
    !	  *			   absorption coefficients.	  *
    !	  * *******************************************************
    !
  
    USE carma_fastjx, ONLY: do3
    USE mem_aerad, ONLY: wave_aerad,u0_aerad, &
  			 sfc_alb_aerad,emisir_aerad,ptop_aerad, &
  			 pbot_aerad,tsfc_aerad
  
    USE mem_globrad, ONLY: nlayer,g, &
  			   ntotal,nprob, &
  			   wave, &
  			   nvert,p,t,o2mol,am, &
  			   co2mol,o3c,o3mol,avg,wol,gol, &
  			   tauray,nsolp,ltemp,nsol,alos,aco2,xaco2, &
  			   ao2,xao2,ao3,xao3,xah2o,psh2o, &
  			   psco2,pso2,pso3,pj,o3mixp, &
  			   akh2o,ako3,akco2,nirp,imie, &
  			   corereal,coreimag
    use chem1_list, only: PhotojMethod
    use mem_chem1,  only: CHEMISTRY
     
    IMPLICIT NONE
   
    INTEGER,INTENT(IN) :: m1,ia,ja,iz,jz
    
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),m1) :: pbar
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),m1) :: o3mix
    INTEGER :: i,i1,j1
    INTEGER :: ii((iz-ia+1)*(jz-ja+1),nlayer)
    INTEGER :: ik((iz-ia+1)*(jz-ja+1),nlayer)
    INTEGER :: ij
    INTEGER :: j
    INTEGER :: k
    INTEGER :: l
    REAL    :: co2mix
    REAL    :: dp((iz-ia+1)*(jz-ja+1),nlayer)
    REAL    :: o2mix
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1))  :: o3mix2
    REAL    :: pm
    REAL    :: ps((iz-ia+1)*(jz-ja+1),nlayer)
    REAL    :: wvo
    REAL    :: x
    INTEGER :: iend
    
    iend=(iz-ia+1)*(jz-ja+1)
  
    !pbar  - layer average pressure (bars)
    !(note - the top layer is from ptop to 0, so average = ptop/2)
    !press - pressure at edge of layer (dyne/cm^2)
    !dpg   - mass of layer (g / cm**2)
    DO ij=1,iend
      pbar(ij,1)  = p_top(ij)/2.0E6
      press(ij,1) = p_top(ij)
    END DO
    DO  k  = 2,nvert
      DO ij=1,iend 
  	pbar(ij,k)  = p_aerad(ij,k-1)/1.0E6
  	press(ij,k) = (p_aerad(ij,k-1) + p_aerad(ij,k)) * 0.5
  	dpg(ij,k-1) = (press(ij,k)-press(ij,k-1)) / g
      END DO
    END DO
    DO ij=1,iend
      pbar(ij,nlayer)  = p_aerad(ij,nvert)/1.0E6
      press(ij,nlayer) = p_surf(ij)
      dpg(ij,nvert)  = (press(ij,nlayer)-press(ij,nvert)) / g
      !skin temperature
      tt(ij,nlayer) = t_surf(ij)
      !amount of water vapor above model domain (gm / cm**2)
      !From 1976 U.S. Standard Atmosphere, mid-latitude sounding:
      !
      !For z_top = 5 km
      !RDH2O(1) = .13
      !
      !For z_top = 2.6 km
      !RDH2O(1) = .64
      rdh2o(ij,1)   = h2ocol_aerad
      !interpolate temperature from layer center (t) to layer edge (tt)
      tt(ij,1) = t_aerad(ij,1)
    END DO
    DO  k = 2, nvert
      DO ij=1,iend 
  	tt(ij,k) = t_aerad(ij,k-1) * (press(ij,k)/p_aerad(ij,k-1)) ** &
  		    (LOG(t_aerad(ij,k)/t_aerad(ij,k-1))/&
  		     LOG(p_aerad(ij,k)/p_aerad(ij,k-1)))

      END DO
    END DO
    !
    !
    !	  DEFINE MASS MIXING RATIOS. O3MIX TAKEN FROM U.S. STANDARD ATMOS-
    !	  PHERE, MID-LATITUDE SOUNDING
    !
    o2mix	   =   0.22*o2mol/am
    co2mix	   =   3.5E-4*co2mol/am
    !
    !	  OZONE COLUMN ABUNDANCE O3C (#/CM**2) ABOVE PTOP WAS CALCULATED
    !	  FROM THE U.S. STANDARD ATMOSPHERE MID-LATITUDE PROFILE.
    !
    !	  This is for z_top = 5 km
    o3c 	   =   9.02E18
    !
    !	  This is for z_top = 3 km
    !	  O3C		 =   9.2E18
    !
    !	  This is for z_top = 1 km
    !	   O3C  	  =   9.3E18
    !
    !	  CONVERT O3C TO MASS MIXING RATIO O3MIX2.
    !
    DO ij=1,iend
       o3mix2(ij) = o3c*o3mol*g/(p_top(ij)*avg)
    END DO
    !  !
    DO  l	    =	nsolp+1,ntotal
      ltemp(l-nsolp)  =   nprob(l) - nsol
    END DO
    !
    x		       =   alos/avg
    !
    !	  CONVERT SOLAR ABSORPTION COEFFICIENTS TO CM**2/GM.
    !
    DO  l	    =	1,nsolp
      !srf	   ACO2(L)	   =   ACO2(L)/(X*CO2MOL)
      !srf	   AO2(L)	 =   AO2(L)/(X*O2MOL)
      !srf	   AO3(L)	 =   AO3(L)/(X*O3MOL)
      aco2(l)	      =   xaco2(l)/(x*co2mol)
      ao2(l)	      =   xao2(l)/(x*o2mol)
      !  if(l.eq.14) print*,'XAO2=',Xao2(l),X,O2MOL
      ao3(l)	      =   xao3(l)/(x*o3mol)
    END DO
    !
    !	  CALCULATE ABSORPTION COEFFICIENTS
    !
    pah2o=0.0
    paco2=0.0
    pao2=0.0
    pao3=0.0
   
    DO  j =   1,nlayer
      DO  l=  1,nsolp
  	DO ij=1,iend
  	  pah2o(ij,l,j)   =  xah2o(l)*pbar(ij,j)**psh2o(l)
  	  paco2(ij,l,j)   =   aco2(l)*pbar(ij,j)**psco2(l)
  	  pao2(ij,l,j)    =   ao2(l)*pbar(ij,j)**pso2(l)
  	  pao3(ij,l,j)    =   ao3(l)*pbar(ij,j)**pso3(l)
  	END DO
      END DO
    END DO
    !
    
    DO ij=1,iend
      DO  j	    =	1,nlayer
  	DO  i	  =   1,6
  	  ii(ij,j)	     =   i
  	  IF(pbar(ij,j) > pj(i)) EXIT
  	END DO
      END DO
    END DO
    
    ps = 0.0
    DO  j	    =	1,nlayer
      DO ij=1,iend
  	IF( ii(ij,j) == 1 ) ii(ij,j) = 2
  	dp(ij,j)	   =   LOG(pj(ii(ij,j)-1)/pj(ii(ij,j)))
  	IF( pbar(ij,j) > pj(6) )THEN
  	  ik(ij,j) = ii(ij,j) - 1
  	ELSE
  	  ik(ij,j) = ii(ij,j)
  	END IF
  	ps(ij,j) = pbar(ij,j)/pj(ik(ij,j))
  	IF (j /= 1) o3mix(ij,j) =o3mixp(ik(ij,j))*ps(ij,j)**(LOG(o3mixp(ii(ij,j)-1)/o3mixp(ii(ij,j)))/dp(ij,j))
      END DO
    END DO
    
    DO  j	=   1,nlayer
      DO  l	  =   1,31
  	DO ij=1,iend
  	  pah2o(ij,nsolp+l,j) = akh2o(l,ik(ij,j))*ps(ij,j)**(LOG (akh2o(l,ii(ij,j)-1)/akh2o(l,ii(ij,j)))/dp(ij,j))
  	END DO
      END DO
      DO  l	  =   32,35
  	DO ij=1,iend
  	  pah2o(ij,nsolp+l,j) = akh2o(32,ik(ij,j))*ps(ij,j)**(LOG (akh2o(32,ii(ij,j)-1)/akh2o(32,ii(ij,j)))/dp(ij,j))
  	  pao3(ij,nsolp+l,j)  = ako3(l-31,ik(ij,j))*ps(ij,j)**(LOG  &
  			     (ako3(l-31,ii(ij,j)-1)/ako3(l-31,ii(ij,j)))/dp(ij,j))
  	END DO
      END DO
      DO ij=1,iend
  	pah2o(ij,nsolp+36,j)  = akh2o(33,ik(ij,j))*ps(ij,j)**(LOG (akh2o(33,ii(ij,j)-1)/akh2o(33,ii(ij,j)))/dp(ij,j))
      END DO
      DO  l	  =   37,40
  	DO ij=1,iend
  	  paco2(ij,nsolp+l,j) = akco2(1,ik(ij,j))*ps(ij,j)**(LOG (akco2(1,ii(ij,j)-1)/akco2(1,ii(ij,j)))/dp(ij,j))
  	  pah2o(ij,nsolp+l,j) = akh2o(l-3,ik(ij,j))*ps(ij,j)**(LOG  &
  		(akh2o(l-3,ii(ij,j)-1)/akh2o(l-3,ii(ij,j)))/dp(ij,j))
  	END DO
      END DO
      DO  l	  =   41,44
  	DO ij=1,iend
  	  paco2(ij,nsolp+l,j)  =  akco2(2,ik(ij,j))*ps(ij,j)**(LOG (akco2(2,ii(ij,j)-1)/akco2(2,ii(ij,j)))/dp(ij,j))
  	  pah2o(ij,nsolp+l,j)  =  akh2o(l-7,ik(ij,j))*ps(ij,j)**(LOG  &
  		(akh2o(l-7,ii(ij,j)-1)/akh2o(l-7,ii(ij,j)))/dp(ij,j))
  	END DO
      END DO
      DO  l	  =   45,48
  	DO ij=1,iend
  	  paco2(ij,nsolp+l,j) = akco2(3,ik(ij,j))*ps(ij,j)**(LOG (akco2(3,ii(ij,j)-1)/akco2(3,ii(ij,j)))/dp(ij,j))
  	  pah2o(ij,nsolp+l,j) = akh2o(l-11,ik(ij,j))*ps(ij,j)**(LOG  &
  	      (akh2o(l-11,ii(ij,j)-1)/akh2o(l-11,ii(ij,j)))/dp(ij,j))
      END DO
    END DO
    DO  l	=   49,51
      DO ij=1,iend
        paco2(ij,nsolp+l,j) = akco2(4,ik(ij,j))*ps(ij,j)**(LOG (akco2(4,ii(ij,j)-1)/akco2(4,ii(ij,j)))/dp(ij,j))
        pah2o(ij,nsolp+l,j) = akh2o(l-11,ik(ij,j))*ps(ij,j)**(LOG  &
  	      (akh2o(l-11,ii(ij,j)-1)/akh2o(l-11,ii(ij,j)))/dp(ij,j))
      END DO
    END DO
    DO  l	=   52,54
      DO ij=1,iend
        paco2(ij,nsolp+l,j) = akco2(5,ik(ij,j))*ps(ij,j)**(LOG (akco2(5,ii(ij,j)-1)/akco2(5,ii(ij,j)))/dp(ij,j))
        pah2o(ij,nsolp+l,j) = akh2o(l-14,ik(ij,j))*ps(ij,j)**(LOG  &
  	      (akh2o(l-14,ii(ij,j)-1)/akh2o(l-14,ii(ij,j)))/dp(ij,j))
      END DO
    END DO
    DO  l	=   55,57
      DO ij=1,iend
        paco2(ij,nsolp+l,j) = akco2(6,ik(ij,j))*ps(ij,j)**(LOG (akco2(6,ii(ij,j)-1)/akco2(6,ii(ij,j)))/dp(ij,j))
        pah2o(ij,nsolp+l,j) = akh2o(l-17,ik(ij,j))*ps(ij,j)**(LOG  &
  	      (akh2o(l-17,ii(ij,j)-1)/akh2o(l-17,ii(ij,j)))/dp(ij,j))
      END DO
    END DO
    DO  l	=   58,nirp
      DO ij=1,iend
  	pah2o(ij,nsolp+l,j) = akh2o(l-17,ik(ij,j))*ps(ij,j)**(LOG  &
  		(akh2o(l-17,ii(ij,j)-1)/akh2o(l-17,ii(ij,j)))/dp(ij,j))
  	END DO
      END DO
    END DO
  
    DO ij=1,iend
      ! store o3mix2 in o3mix(1)
      o3mix(ij,1) = o3mix2(ij)
    END DO
    
    !	  here we find taugas. it is tauco2+tauo2+tauo3.
    DO  l=1,ntotal
      DO ij=1,iend
  	pm=p_top(ij)/g
  	taugas(ij,l,1) = pm*(o2mix*pao2(ij,l,1)+co2mix* &
  			     paco2(ij,l,1)+o3mix(ij,1)*pao3(ij,l,1)) 
      END DO
    END DO
    DO    j =	2,nlayer
      DO  l    =   1,ntotal
  	 DO ij=1,iend
  	   pm=dpg(ij,j-1)
  	   taugas(ij,l,j) = pm*(o2mix*pao2(ij,l,j)+co2mix* &
  			      paco2(ij,l,j)+o3mix(ij,j)*pao3(ij,l,j)) 
  	 END DO
      END DO
    END DO
  
    !
    !	  wave must be in microns
    !	  calculate rayleigh optical depth parameters.
    !
    DO  l      =   1,ntotal
      wvo	=    wave(nprob(l))
      tauray(l) =   (8.46E-9/wvo**4) * ( 1.+0.0113/wvo**2+0.00013/wvo**4 )
    END DO
  
    !	  we do not include rayleigh scattering in infrared
    DO  j = 1,nvert
      DO  l= 1,ntotal
  	DO ij=1,iend
  	  IF( l <= nsolp ) THEN
  	    paray(ij,l,j+1) = tauray(l)*dpg(ij,j)*g
  	  ELSE
  	    paray(ij,l,j+1) = 0.
  	  END IF
  	END DO
      END DO
      DO  l   =   1,ntotal
  	DO ij=1,iend
  	  IF( l <= nsolp ) THEN
  	    paray(ij,l,1) = tauray(l)*p_top(ij)
  	  ELSE
  	    paray(ij,l,1) = 0.
  	  END IF
  	END DO
      END DO
    END DO

!-srf- 06-05-2008 - using Harrington data set
!- unit of "do3" is g/g, needs to be changed to molec/cm^2 at fastjx routine
!    if (trim(PhotojMethod) == 'FAST-JX' .and. CHEMISTRY > 0) & 
!    CALL C_2d_3d(o3mix,do3(:,ia:iz,ja:jz),m1,ia,iz,ja,jz,iend) !LFR
   
  END SUBROUTINE setuprad
  
  SUBROUTINE calcproperties
    ! **********************************************************************
    !
    !		 CALCULATE THE AEROSOL EXTINCTION CROSS SECTIONS
    !
    ! **********************************************************************
    !
    !	  Get <is_grp_ice> and radius grid from interface common block
    !	  and calculate cross-sectional area for each bin.
    !
    USE mem_aerad, ONLY: is_grp_ice_aerad,r_aerad,rcore_aerad,rup_aerad, &
  			 rcoreup_aerad,lunmie,lunoprt,ir_above_aerad, &
  			 tabove_aerad
  
    USE mem_globrad, ONLY: ngroup,nrad,core_rad,xsecta,pi, &
  			   coreup_rad,i_write,i_read,nwave,rmin,rdqext, &
  			   qscat,qbrqs,nsol,wave,corereal,coreimag,nsolp, &
  			   weight,ntotal,sol,solfx,nprob,iblackbody_above, &
  			   t_above,ncount,nlow,plank,sbk
  
    USE mem_globaer, ONLY: r,rcore,rup,rcoreup,is_grp_ice
  
    IMPLICIT NONE
    
    INTEGER :: i
    INTEGER :: ibeyond_spectrum
    INTEGER :: ig
    INTEGER :: i_mie
    INTEGER :: irefr
    INTEGER :: j
    INTEGER :: jj
    INTEGER :: k
    INTEGER :: l
    INTEGER :: mgroup
    INTEGER :: mrad
    INTEGER :: mwave
    INTEGER :: n_thetd
    LOGICAL :: all_ok
    REAL    :: awave
    REAL    :: corerad
    REAL    :: ctbrqs
    REAL    :: ddr
    REAL    :: ddrc
    REAL    :: qextd
    REAL    :: qscatd
    REAL    :: r_real
    REAL    :: rr
    REAL    :: sum
    REAL    :: sum1
    REAL    :: sum2
    REAL    :: t1
    REAL    :: thetd(1)
    REAL    :: tmag
    REAL    :: v
    REAL    :: wvno
  
    CHARACTER(LEN=*),PARAMETER :: &
    lab355='(//,"setuprad: error in weights ",/," ' &
  	   //'sum of weights for solar =",1pe15.5,/,"' &
  	   //' sum of weights for ir = ",1pe15.5,/," ' &
  	   //'total sum =  ",1pe15.5)'
    
    DO ig = 1, ngroup
      DO I = 1, nrad
  	xsecta(i,ig) = pi * r(i,ig)**2.
      END DO
    END DO
    !
    !	  Set <i_mie> = I_READ to WRITE the mie coefficients to a data file,
    !		      = I_WRITE to READ them
    !
  	   i_mie = I_READ
    !i_mie = i_WRITE
    !srf - nao precisa escrever o arquivo mie.data
    !	   i_mie = 555
   
    IF ( i_mie .eq. i_READ ) THEN
      !   READ extinction and scattering coefficients
      OPEN(lunmie,file='./mie.data',form='FORMATted')
      !   Check that input file is consistent with radius and
      !   wavelength grids
      !
      READ(lunmie,*) mwave,mrad,mgroup
      all_ok = (mwave .eq. nwave) .and. (mrad .eq. nrad) .and. (mgroup .eq. ngroup)
      !
      READ(lunmie,*) (rmin(i),i=1,mgroup)
      !
      DO ig = 1, ngroup
  	!...dbg:
  	all_ok = all_ok .and. ( r(1,ig) .eq. rmin(ig) )
      END DO
      !
      IF ( .not. all_ok )THEN
  	WRITE(lunoprt,*) ' setuprad: mie.data grid(s) bad: '
  	WRITE(lunoprt,*) ' in mie.data, mwave, mrad, mgroup = ',mwave, mrad, mgroup
  	WRITE(lunoprt,*) ' and rmin = ', rmin
  	stop 1
      END IF
      DO ig = 1,ngroup
  	DO i = 1,nrad
  	  DO l = 1,nwave
  	    READ(lunmie,*) rdqext(i,ig,l), qscat(i,ig,l), qbrqs(i,ig,l)
  	  END DO
  	END DO
      END DO
      CLOSE(lunmie)
    ELSE
      ! calculate extinction and scattering coefficients
      !
      DO ig = 1,ngroup
  	!
  	!	select ice/liquid index of refractive index array
  	!
  	IF( is_grp_ice(ig) )THEN
  	  irefr = 2
  	ELSE
  	  irefr = 1
  	END IF
	!
  	!	<thetd> is angle between incident and scattered radiation
  	!	<j_thetd> is number of <thetd> values to consider
  	!
  	thetd = 0.0
  	n_thetd = 1
  	DO  l=1,nwave
  	  !
  	  !kml for biomass burning particles:
  	  r_real = 1.495
  	  tmag = 1.0e-6
  
  	  !	  real = 1.52
  	  !	  tmag = 0.015

  	  !
  	  !	 calculate the center of the wavelength interval of an ir interval
  	  !
  	  IF( l .le. nsol ) THEN
  	    awave = wave(l)
  	  ELSE
  	    awave = 0.5*(wave(l)+wave(l+1))
  	  END IF
  	  wvno  	    =	2.*pi/(awave*1.0e-4)
  	  !
  	  DO i=1,nrad
  	    IF(i .eq. 1) THEN
  	      ddr		 =   0.2*(rup(1,ig)-r(1,ig))
  	      rr		 =   r(1,ig)
  	      corerad	 =   rcore(1,ig)
  	      ddrc	 =   0.2*(rcoreup(i,ig)- rcore(1,ig))
  	    ELSE
  	      ddr		 =   0.2*(rup(i,ig)-rup(i-1,ig))
  	      rr		 =   rup(i-1,ig)
  	      corerad	 =   rcoreup(i-1,ig)
  	      ddrc	 =   0.2*(rcoreup(i,ig)-rcoreup(i-1,ig))
  	    END IF
  	    !
  	    rdqext(i,ig,l)   =   0.0
  	    qscat(i,ig,l)    =   0.0
  	    qbrqs(i,ig,l)    =   0.0
  	    !
  	    DO j=1,6
  	      !
  	      !
  	      ! limit x=2*pi/wave to no larger 1000 to avoid anguish in the mie code.
  	      !
  	      IF( wvno*rr .gt. 1000. ) rr = 1000./wvno
  	      !  print*,'-------------------------------------------------'
  	      !     IF(j.eq.1) rr = 0.1 * 1.e-4
  	      !    corerad=0.5*rr
  	      !  print*, rr,corerad,rr/corerad
  	      !  print*, rr,real,tmag,thetd,n_thetd,qextd,qscatd,ctbrqs,
  	      ! 1	       corerad,corereal,coreimag,wvno
  	      !  print*,'-------------------------------------------------'
  	      !      stop
  	      CALL miess(rr,r_real,tmag,thetd,n_thetd,qextd,qscatd,ctbrqs,corerad,corereal,coreimag,wvno)
  	      !     IF(l.eq.7 .or. l.eq.8)
  	      ! &      print*,'qex qsc=',qextd,qscatd
  	      !
  	      rdqext(i,ig,l)	 =   rdqext(i,ig,l)+qextd/6.
  	      qscat(i,ig,l)	 =   qscat(i,ig,l)+qscatd/6.
  	      qbrqs(i,ig,l)	 =   qbrqs(i,ig,l)+ctbrqs/6.
  	      rr		     =   rr+ddr
  	      corerad	     =   corerad + ddrc
  	      !
  	    END DO
  	  END DO
  	END DO
  	!
  	!	 stop
      END DO	  ! ig=1,ngroup
      !
    END IF
  
    !
    !srf - nao precisa escrever o arquivo mie.data
    !	 i_mie=555
  
    !	   DO ig = 1,ngroup
    !	    DO i = 1,nrad
    !	     DO l = 1,nwave
    !	       print*, 'rdqext(i,ig,l), qscat(i,ig,l), qbrqs(i,ig,l)'
    !	       print*, rdqext(i,ig,l), qscat(i,ig,l), qbrqs(i,ig,l)
    !	     END DO
    !	    END DO
    !	   END DO
  
    !IF ( i_mie .eq. i_WRITE ) THEN
    !  PRINT *,'LFR->006.1.4: Writing mie';CALL flush(6)
   !
   !   ! WRITE extinction and scattering coefficients to data file
   !   !
   !   OPEN(lunmie,file='./mie.data',form='formatted')
   !   !
   !   PRINT *,'LFR->006.1.4.1: Writing nwave,nrad,ngroup',nwave,nrad,ngroup;CALL flush(6)
   !   WRITE(lunmie,*) nwave,nrad,ngroup
   !   PRINT *,'LFR->006.1.4.1: Writing r';CALL flush(6)
   !   DO ig = 1,ngroup
   !	 DO i = 1,nrad
  !	  WRITE(lunmie,*) r(i,ig)
  !	 END DO
  !    END DO
  !
  !    PRINT *,'LFR->006.1.4.2: Writing rdqext,qscat,qbrqs';CALL flush(6)
  !    DO ig = 1,ngroup
  !	 DO i = 1,nrad
  !	  DO l = 1,nwave
  !	    WRITE(lunmie,*) rdqext(i,ig,l), qscat(i,ig,l), qbrqs(i,ig,l)
  !	  END DO
  !	 END DO
  !    END DO
  !
  !    CLOSE(lunmie)
  !    PRINT *,'LFR->006.1.4.3: Mie close';CALL flush(6)
  !
  !  END IF
    !	 stop
    !
    !	WRITE some values to print file
    !
    !	 DO ig = 1, ngroup
    !	   WRITE (lunoprt,500) ig
    !	   DO i = 1, nrad
    !	     sizparm6 = 2.*pi*r(i,ig)/(wave(6)*1.e-4)
    !	     sizparm24 = 2.*pi*r(i,ig)/(wave(24)*1.e-4)
    !	     WRITE(lunoprt,505) i,r(i,ig),rdqext(i,ig,6),sizparm6,
    !	1	rdqext(i,ig,24),sizparm24
    !	   END DO
    !	 END DO
    !
    !500  FORMAT(/," setuprad: igroup = ",i4,//,"   i	 r(cm)   rdqext(6)	x6","	rdqext(24)	x24",/)
    !505  FORMAT(i4,5(1pe11.2))
    !
    ! *********************************************************************
    !
    !				 check sum of weights
    !
    ! **********************************************************************
    !
    sum 	=   0.0
    sum1	=   0.0
    sum2	=   0.0
    DO l	=   1,nsolp
      sum	     =   sum+weight(l)
    END DO
    DO l	=   nsolp+1,ntotal
      sum1	 =   sum1+weight(l)
    END DO
    sum2	=   sum+sum1
  
    !
    IF ( abs(nwave-sum2) .gt. 1.e-3 ) WRITE(lunoprt,FMT=lab355) sum,sum1,sum2
    !
    DO  l   =	1,nsolp
      sol(l)  =   solfx(nprob(l)) * weight(l)
    END DO
    !	 print*, 'wave(l),nprob(l),weight(l),solfx(nprob(l)),sol(l)'
    !	 DO 361 l   =	1,ntotal
    !
    !	   print*,wave(nprob(l)),nprob(l),weight(nprob(l)),solfx(nprob(l)),sol(l)
    !
    ! 361  continue
  
    !
    ! *********************************************************************
    !
    !	compute planck function table. wave is in units of microns.
    !
    ! **********************************************************************
    !
    !	set <iblackbody_above> = 1 to include a source of radiation
    !	at the top of the radiative transfer model DOmain
    
    iblackbody_above = ir_above_aerad
    t_above = tabove_aerad
    !
    !	set <ibeyond_spectrum> = 1 to include blackbody radiation at
    !	wavelengths longer than wave(nwave+1) in plank(nwave+1-nsol) and
    !	at wavelengths shorter than wave(nsol+1) in plank(1)
    !
    ibeyond_spectrum = 1
    !
    IF( ibeyond_spectrum .eq. 1 )THEN 
      DO j  =	1,ncount
  	plank(nwave+1-nsol,j) = (0.01*float(nlow+j))**4
      END DO
      DO i =   nsol+2,nwave
        DO j  =	1,ncount
  	  k =	i-nsol
  	  v =	1.438e4  /  wave(i)
  	  CALL plnk(v,(0.01*float(nlow+j)),plank(k,j))
  	END DO
      END DO
    ELSE 
      DO i =   nsol+1,nwave+1
        DO j  =	1,ncount
  	  k =	i-nsol
  	  v =	1.438e4  /  wave(i)
  	  CALL plnk(v,(0.01*float(nlow+j)),plank(k,j))
  	END DO
      END DO
    END IF
    !
    DO j   =   1,ncount
  
      IF( ibeyond_spectrum .eq. 1 )THEN
  
  	plank(1,j) = plank(2,j)*sbk/pi
  	DO l  =   nsol+2,nwave
  	  k  =   l-nsol
  	  plank(k,j) = (plank(k+1,j)-plank(k,j))*sbk/pi
  	END DO
  
      ELSE
  
  	DO l  =   nsol+1,nwave
  	  k  =   l-nsol
  	  plank(k,j) = (plank(k+1,j)-plank(k,j))*sbk/pi
  	END DO
  
      END IF
  
    END DO
    !
  
  END SUBROUTINE calcproperties
  
!KML2!!!!!!!!1
  SUBROUTINE nocalcproperties
    ! **********************************************************************
    !
    !		 CALCULATE THE AEROSOL EXTINCTION CROSS SECTIONS
    !
    ! **********************************************************************
    !
    !	  Get <is_grp_ice> and radius grid from interface common block
    !	  and calculate cross-sectional area for each bin.
    !
    USE mem_aerad, ONLY: ir_above_aerad,tabove_aerad,lunoprt


  
    USE mem_globrad, ONLY: ngroup,nrad,xsecta,pi, &
  			   i_write,i_read,nwave,rmin, &
  			   nsol,wave,nsolp, &
  			   weight,ntotal,sol,solfx,nprob,iblackbody_above, &
  			   t_above,ncount,nlow,plank,sbk
  
    USE mem_globaer, ONLY: r
  
    IMPLICIT NONE
    
    INTEGER :: i
    INTEGER :: ibeyond_spectrum
    INTEGER :: ig
    INTEGER :: j
    INTEGER :: jj
    INTEGER :: k
    INTEGER :: l
    REAL    :: sum
    REAL    :: sum1
    REAL    :: sum2
    REAL    :: v
  
    CHARACTER(LEN=*),PARAMETER :: &
    lab355='(//,"setuprad: error in weights ",/," ' &
  	   //'sum of weights for solar =",1pe15.5,/,"' &
  	   //' sum of weights for ir = ",1pe15.5,/," ' &
  	   //'total sum =  ",1pe15.5)'
    
    DO ig = 1, ngroup
      DO I = 1, nrad
  	xsecta(i,ig) = pi * r(i,ig)**2.
      END DO
    END DO
   

    !
    ! *********************************************************************
    !
    !				 check sum of weights
    !
    ! **********************************************************************
    !
    sum 	=   0.0
    sum1	=   0.0
    sum2	=   0.0
    DO l	=   1,nsolp
      sum	     =   sum+weight(l)
    END DO
    DO l	=   nsolp+1,ntotal
      sum1	 =   sum1+weight(l)
    END DO
    sum2	=   sum+sum1
  
    !
    IF ( abs(nwave-sum2) .gt. 1.e-3 ) WRITE(lunoprt,FMT=lab355) sum,sum1,sum2
    !
    DO  l   =	1,nsolp
      sol(l)  =   solfx(nprob(l)) * weight(l)
    END DO


    ! *********************************************************************
    !
    !	compute planck function table. wave is in units of microns.
    !
    ! **********************************************************************
    !
    !	set <iblackbody_above> = 1 to include a source of radiation
    !	at the top of the radiative transfer model DOmain
    
    iblackbody_above = ir_above_aerad
    t_above = tabove_aerad
    !
    !	set <ibeyond_spectrum> = 1 to include blackbody radiation at
    !	wavelengths longer than wave(nwave+1) in plank(nwave+1-nsol) and
    !	at wavelengths shorter than wave(nsol+1) in plank(1)
    !
    ibeyond_spectrum = 1
    !
    IF( ibeyond_spectrum .eq. 1 )THEN 
      DO j  =	1,ncount
  	plank(nwave+1-nsol,j) = (0.01*float(nlow+j))**4
      END DO
      DO i =   nsol+2,nwave
        DO j  =	1,ncount
  	  k =	i-nsol
  	  v =	1.438e4  /  wave(i)
  	  CALL plnk(v,(0.01*float(nlow+j)),plank(k,j))
  	END DO
      END DO
    ELSE 
      DO i =   nsol+1,nwave+1
        DO j  =	1,ncount
  	  k =	i-nsol
  	  v =	1.438e4  /  wave(i)
  	  CALL plnk(v,(0.01*float(nlow+j)),plank(k,j))
  	END DO
      END DO
    END IF
    !
    DO j   =   1,ncount
  
      IF( ibeyond_spectrum .eq. 1 )THEN
  
  	plank(1,j) = plank(2,j)*sbk/pi
  	DO l  =   nsol+2,nwave
  	  k  =   l-nsol
  	  plank(k,j) = (plank(k+1,j)-plank(k,j))*sbk/pi
  	END DO
  
      ELSE
  
  	DO l  =   nsol+1,nwave
  	  k  =   l-nsol
  	  plank(k,j) = (plank(k+1,j)-plank(k,j))*sbk/pi
  	END DO
  
      END IF
  
    END DO
    !
  
  END SUBROUTINE nocalcproperties


!KML2!!!!!!!!

  SUBROUTINE prerad(m1,dztr,fmapt,ia,iz,ja,jz,nzpmax,m2,m3)
    !
    !  Note that vertical index in radiative transfer model DOmain is reversed
    !  for cartesian coordinates.
    !
    !  Indices <ix> and <iy> are passed through global COMMON block.
    !
    !  INCLUDE global constants and variables
    !
    
    !USE mem_aerad, ONLY: qv_aerad,pc_aerad
    USE mem_aerad, ONLY: nbin
    USE mem_globaer,ONLY: ngroup,ienconc,nelem
    
    IMPLICIT NONE
  
    INTEGER,INTENT(IN)  	      :: m1,m2,m3,ia,iz,ja,jz,nzpmax
    REAL   ,INTENT(IN), DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax) :: dztr
    REAL   ,INTENT(IN), DIMENSION((iz-ia+1)*(jz-ja+1)) :: fmapt
  
    INTEGER :: ibin
    INTEGER :: i,j,ij,iend
    INTEGER :: iep
    INTEGER :: igas
    INTEGER :: igroup
    INTEGER :: k,kk
    INTEGER :: nzz
    REAL    :: xymet
  
    iend=(iz-ia+1)*(jz-ja+1)
  	   
    !  Load profiles of temperature [K], water vapor [g/g], and
    !  aerosol particle concentrations [#/cm^2]
    !
    igas = 1
    !srf
    nzz = m1 - 1
  
    DO k = 1,NZZ
       !  Reverse the vertical index when in cartesian coordinates
       kk = nzz + 1 - k
       ! For radiation code: qv-aerad have g[H20]/g[ar] unit
       DO ij=1,iend
  	 qv_aerad(ij,kk) = gc(ij,k,igas) / rhoa(ij,k)
       END DO
       DO igroup = 1,ngroup
  	  iep = ienconc(igroup)
  	  DO ibin = 1,nbin
  	    DO ij=1,iend 
  	      xymet = fmapt(ij)*fmapt(ij)
  	      pc_aerad(ij,kk,ibin,igroup) = pc(ij,k,ibin,iep) *  &
  		  (1./dztr(ij,k)) / xymet
  	    END DO
  	  END DO
  	END DO
    END DO
    !
  
  END SUBROUTINE prerad
!kmlnew  
  SUBROUTINE radtran(albedt,cosz,m1,m2,m3,ia,iz,ja,jz,aot11)
!kmlnew     
    USE mem_aerad, ONLY: u0_aerad,qrad_aerad, &
  			 alb_toai_aerad,alb_tomi_aerad,alb_toa_aerad, &
  			 fsl_up_aerad,fsl_dn_aerad,fir_up_aerad,fir_dn_aerad, &
                         nir
      
    USE mem_globrad, ONLY: isl,nvert,nlayer, &
  			   ngroup,nrad,u0,nsolp,albedo_sfc, &
  			   emis,ntotal,emisir,ir,irs, &
  			   fdegday,g,scday,qrad,pi,epsilon,xsecta,rdqext, &
  			   nprob,qscat,nwave,tslu,tsld,fupbs,fdownbs, &
  			   fnetbs,nsol,fslu,fsld,alb_toa,alb_tomi,alb_toai, &
  			   solfx, &
  			   tiru,fupbi,fdownbi,fnetbi,firu,xirup
  
    IMPLICIT NONE
  
    INTEGER,INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
    REAL,INTENT(IN),DIMENSION((iz-ia+1)*(jz-ja+1)) :: albedt,cosz
    REAL  :: aot11((iz-ia+1)*(jz-ja+1))
    INTEGER :: i,i1,j1
    INTEGER :: ig
    INTEGER :: j
    INTEGER :: l,k
    REAL    :: term1((iz-ia+1)*(jz-ja+1))
    INTEGER :: count=0
    INTEGER :: ij,iend
    REAL :: heati((iz-ia+1)*(jz-ja+1),nlayer)
    REAL :: heats((iz-ia+1)*(jz-ja+1),nlayer)
    REAL :: heat((iz-ia+1)*(jz-ja+1),nlayer)
    
    iend=(iz-ia+1)*(jz-ja+1)
  
    heats   =  0.0
    heati   =  0.0
    
    !
    !	  interpolate temperature from layer center (t) to layer edge (tt)
    DO ij=1,iend
      tt(ij,1) = t_aerad(ij,1)
    END DO
    DO  j = 2, nvert
      DO ij=1,iend
  	  tt(ij,j) = t_aerad(ij,j-1) * (press(ij,j)/p_aerad(ij,j-1)) ** &
  			(LOG(t_aerad(ij,j)/t_aerad(ij,j-1))/ &
  		    LOG(p_aerad(ij,j)/p_aerad(ij,j-1)))
      END DO
    END DO
  
    !	  water vapor (g / cm**2)
      DO  j = 2, nlayer
      DO ij=1,iend
  	  rdh2o(ij,j)	= qv_aerad(ij,j-1) * dpg(ij,j-1)
      END DO
    END DO
  
    !	  aerosol concentrations (# / cm**2)
    DO ig = 1, ngroup
      DO  j = 2, nvert
  	DO  i = 1, nrad
  	  DO ij=1,iend
  	    caer(ij,j,i,ig)  = pc_aerad(ij,j-1,i,ig)  !!!!!!
  	  END DO
  	END DO
      END DO
    END DO  
  
    !surface reflectivity and emissivity
    DO  l =  1,nsolp
      DO ij=1,iend
  	rsfx(ij,l) =  albedt(ij)
  	emis(l) =  0.0
      END DO
    END DO
    DO  l =  nsolp+1,ntotal
      DO ij=1,iend
  	emis(l) =  emisir_aerad
  	rsfx(ij,l) = 1.0 - emis(l)
      END DO
    END DO
    
    !set wavelength limits lla and lls based on values of isl and ir
    lla=  ntotal
    lls=  1
  
    DO ij=1,iend
      IF(isl_aerad(ij)  == 0) THEN
  	lls(ij)   =  nsolp+1
      END IF
    END DO
    !
    IF(ir_aerad   == 0) THEN
      DO ij=1,iend
  	lla(ij)  =  nsolp
      END DO
    END IF
     
      !DO ij=1,iend
      !     print*,'AOT11 na radtran=',ij,aot11
      !END DO
	   
    !calculate the optical properties
    CALL oppr(ia,iz,ja,jz,m1,aot11)
    !
    !	  if infrared calculations are required then calculate
    !	  the plank function
    !
    IF(ir_aerad /= 0) THEN
      CALL oppr1(ia,iz,ja,jz,m1)
    END IF
    !
    !	  if no infrared scattering then set index to number of
    !	  solar intervals
    !
    IF(irs == 0) THEN
      lla  =  nsolp
    END IF
    !
    !	  if either solar or infrared scattering calculations are required
    !	  call the two stream code and find the solution
    !
    
  
    CALL twostr(m1,ia,iz,ja,jz)
    
    !DO i1=ia,iz
    !  DO j1=ja,jz
    !	 IF(isl_aerad(i1,j1) /= 0 .OR. irs .NE. 0 ) THEN
    !	   CALL add(i1,j1,ia,iz,ja,jz,cosz,m2,m3)
    !	 END IF
    !  END DO
    !END DO
    CALL add(m1,ia,iz,ja,jz,cosz,m2,m3)
  
    !
    !	  if infrared calculations are required then call newflux1 for
    !	  a more accurate solution
    !
    IF(ir_aerad /= 0) THEN
      CALL newflux1(m1,ia,iz,ja,jz)
    END IF
    
    !	  calculate infrafred and solar heating rates (deg/day),
    DO  j      =  1,nvert
      DO ij=1,iend
  	IF(isl_aerad(ij) /= 0) THEN
  	  term1(ij)	 =  fdegday/(dpg(ij,j)*g)
  	END IF
      END DO
      DO  l =  1,nsolp
  	DO ij=1,iend
  	  IF(isl_aerad(ij) /= 0) THEN
  	      heats(ij,j)   =  heats(ij,j)+(fnet(ij,l,j+1)-fnet(ij,l,j))*term1(ij)
  	  END IF
  	END DO
      END DO
    END DO
    !  
    DO  j      =  1,nvert
      DO ij=1,iend
  	IF(ir_aerad /= 0) THEN
  	  term1(ij)	 =  fdegday/(dpg(ij,j)*g)
  	END IF
      END DO
      DO  l =  nsolp+1,ntotal
  	DO ij=1,iend
  	  IF(ir_aerad /= 0) THEN
  	    heati(ij,j)  =  heati(ij,j)+(directu(ij,l,j+1)-direc(ij,l,j+1)  &
  		       -(directu(ij,l,j)-direc(ij,l,j)) )*term1(ij)
	  
!srf
!            if(heati(ij,j) < -10) then
!	      print*,ij,j,heati(ij,j)
!        if((ii == 41 .and. jj==36) .or.(ii == 42 .and. jj==36) ) then
!	      if(l==nsolp+1)print*,'radtran ii jj l j',ii,jj,l,j
!	      print*,j,directu(ij,l,j+1),direc(ij,l,j+1)  &
 ! 		       ,directu(ij,l,j),direc(ij,l,j),term1(ij)
!	     endif
!srf		       
  
	  
	  
  	  END IF
  	END DO
      END DO
    END DO
    !
    DO j      =  1,nvert
      DO ij=1,iend
  	!     Load heating rates [deg_K/s] into interface common block
  	heat(ij,j)	   =  heats(ij,j)+heati(ij,j)
  	heats_aerad(ij,j) =  heats(ij,j)/scday
  	heati_aerad(ij,j) =  heati(ij,j)/scday
      END DO
    END DO
    
    !DO  j	=  1,nvert
    ! heats(j)   =  0.0
    !  term1	  =  fdegday/(dpg(j)*g)
    !
    !  IF(isl /= 0) THEN
    !	 DO  l    =  1,nsolp
    !	  heats(j)   =  heats(j)+(fnet(l,j+1)-fnet(l,j))*term1
    !		    print*,l,heats(j)
    !	 END DO
    !  END IF
    !
    !  IF(ir /= 0) THEN
    !	 DO  l    =  nsolp+1,ntotal
    !	  heati(j)  =  heati(j)+( directu(l,j+1)-direc(l,j+1)  &
    !	      -(directu(l,j)-direc(l,j)) )*term1
    !	 END DO
    !  END IF
    !  heat(j)       =  heats(j)+heati(j)
    !
    !	  Load heating rates [deg_K/s] into interface common block
    !
    !  heats_aerad(j) =  heats(j)/scday
    !  heati_aerad(j) =  heati(j)/scday
    !
    !END DO
    !
    !	  Here we Calculate (4 * pi * mean_intensity) for the IR.
    !
    !DO j = 1, nvert
    !  DO l = nsolp+1, ntotal
  !	DO ij=1,iend
  !	  IF (ir_aerad /= 0) THEN
  !	    tmi(ij,l,j) = tmiu(ij,l,j)+tmid(ij,l,j)
  !	  END IF
  !	END DO
  !    END DO
  !  END DO
    !
    !	  Here we compute the heating rates for droplets
    !	  (C11 converts W/m^2 to erg/cm^2)
    !
    !qrad = 0.
    !c11 = 1000.
    !IF (ir /= 0) THEN
    !  DO ig = 1, ngroup
    !	 DO i = 1, nrad
    !	   DO j = 1, nvert
    !	    DO l = nsolp+1, ntotal
    !	      x = tmi(l,j)-4.0*pi*ptemp(l,j)
    !	      IF( ABS(x/tmi(l,j)) < epsilon ) x = 0.
    !	      qrad(i,j,ig) = qrad(i,j,ig) + x*c11*xsecta(i,ig) *  &
    !		  (rdqext(i,ig,nprob(l))-qscat(i,ig,nprob(l)))
    !	    END DO
    !	   END DO  ! j=1,nvert
    !	 END DO   ! i=1,nrad
    !  END DO	    ! ig=1,ngroup
    !END IF
    !
    !IF (isl /= 0) THEN
    !  DO ig = 1, ngroup
    !	 DO i = 1, nrad
    !	   DO j = 1, nvert
    !	    DO l = 1, nsolp
    !	      qrad(i,j,ig) = qrad(i,j,ig) + tmi(l,j)*c11*xsecta(i,ig) *  &
    !		  (rdqext(i,ig,nprob(l))-qscat(i,ig,nprob(l)))
    !	    END DO
    !	   END DO  ! j=1,nvert
    !	 END DO   ! i=1,nrad
    !  END DO	    ! ig=1,ngroup
    !END IF
    !
    !	  Load layer averages of droplet heating rates into interface common block
    !
    !DO ig = 1, ngroup
    !  DO i = 1, nrad
    !	 DO j = 1, nvert
    !	   IF (j == nvert) THEN
    !	    qrad_aerad(i,j,ig) = qrad(i,j,ig)
    !	  ELSE IF (j > 1) THEN
    !	    qrad_aerad(i,j-1,ig) = 0.5 * ( qrad(i,j,ig) + qrad(i,j-1,ig) )
    !	  END IF
    !	 END DO    ! j=1,nvert
    !  END DO	  ! i=1,nrad
    !END DO	  ! ig=1,ngroup
    !
    !
    !	Calculate some diagnostic quantities (formerly done in radout.f) and
    !	load them into the interface common block.  None of these presently
    !	influence any microphysical processes -- hence, the following code only
    !	needs to be executed before the aerosol model writes its output.
    !	Not all of the calculated quantities are presently being
    !	loaded into the interface common block.
    !
    !	Load optical depths into interface common block
    !
    !DO i = 1, nwave
    !  opdaerad(i) = uopd(i,nlayer)
    !	print*,'AOT',uopd(i,nlayer),opd(i,nlayer)
    !END DO
   
    !
    !	  <tsLu> and <tsLd> are total upwelling and downwelling solar
    !	  fluxes at top-of-atmosphere
    !
    !tslu = 0.
    !tsld = 0.
    !
    !	  <fupbs>, <fdownbs>, and <fnetbs> are total upwelling, downwelling, and net
    !	  solar fluxes at grid boundaries
    !
    !DO  j = 1, nlayer
    !  fupbs(j) = 0.
    !  fdownbs(j) = 0.
    !  fnetbs(j) = 0.
    !END DO
    !
    !	  <fsLu> and <fsLd> are upwelling, downwelling, and net
    !	  solar fluxes at top-of-atmosphere (spectrally-resolved)
    !
    !	  <alb_toa> is albedo at top-of-atmosphere (spectrally-resolved)
    !
    !DO  i = 1, nsol
    !  fslu(i) = 0.
    !  fsld(i) = 0.
    !  alb_toa(i) = 0.
    !END DO
    !
    !	   <alb_tomi> and <alb_toai> are total solar albedos at top-of-model
    !	   and top-of-atmosphere
    !
    !alb_tomi = 0.
    !alb_toai = 0.
    !
    !	  CALCULATE SOLAR ABSORBED BY GROUND, SOLNET, AND UPWARD AND DOWNWARD
    !	  LONGWAVE FLUXES AT SURFACE, XIRUP AND XIRDOWN (WATTS/M**2)
    !
    solnet  = 0.0
   
    DO  l   =  1,nsolp
      DO ij=1,iend
  	IF (isl_aerad(ij) /= 0) THEN
  	  solnet(ij) = solnet(ij) - fnet(ij,l,nlayer)
    !fp = ck1(l,1)*el2(l,1) - ck2(l,1)*em2(l,1) + cp(l,1)
    !fslu( nprob(l) ) = fslu( nprob(l) ) + fp
    !DO  j = 1, nlayer
    !	  fp = ck1(l,j)*el1(l,j) + ck2(l,j)*em1(l,j) + cpb(l,j)
    !	  fupbs(j) = fupbs(j) + fp
    !	  fnetbs(j) = fnetbs(j) + fnet(l,j)
    !	  IF (l == nsolp) fdownbs(j) = fupbs(j) - fnetbs(j)
    !	 END DO
  	END IF
      END DO
    END DO
      !DO  i = 1, nsol
  	!fsld(i) = u0*solfx(i)
  	!alb_toa(i) = fslu(i)/fsld(i)
  	!tslu = tslu + fslu(i)
  	!tsld = tsld + fsld(i)
      !END DO
    !
      !alb_tomi = fupbs(1)/fdownbs(1)
      !alb_toai = tslu/tsld
    !
    !	   Load albedos into interface common block
    !
      !alb_toai_aerad = alb_toai
      !alb_tomi_aerad = alb_tomi
      !DO i = 1, nsol
      !  alb_toa_aerad(i) = alb_toa(i)
      !END DO
    !
    !	   Load fluxes into interface common block
    !
      !DO j = 1, nlayer
      !  fsl_up_aerad(j) = fupbs(j)
      !  fsl_dn_aerad(j) = fdownbs(j)
      !END DO
   
   ! END DO
    !
    !	  <tiru> is total upwelling infrared flux at top-of-atmosphere;
    !	  <fupbi>, <fdownbi>, and <fnetbi> are total upwelling, downwelling, and net
    !	  infrared fluxes at grid boundaries
    !
    !tiru = 0.
    !DO  j = 1, nlayer
    !  fupbi(j)   =  0.0
    !  fdownbi(j)   =  0.0
    !  fnetbi(j)  =  0.0
    !END DO
    !
    !	  <firu> is upwelling infrared flux at top-of-atmosphere (spectrally-resolved)
    !
    !DO  i = 1, nir
    !  firu(i) = 0.
    !END DO
   
    xirdown = 0.0
    !xirup   = 0.0
   
    DO  l  =  nsolp+1,ntotal
      DO ij=1,iend
  	IF (ir_aerad /= 0) THEN
  	    xirdown(ij) = xirdown(ij) + direc(ij,l,nlayer)

    !	 xirup   = xirup  + directu(l,nlayer)
    !	 firu( nprob(l)-nsol ) = firu( nprob(l)-nsol ) + directu(l,1)
    !	 DO  j = 1, nlayer
    !	  fupbi(j) = fupbi(j) + directu(l,j)
    !	  fdownbi(j) = fdownbi(j) + direc  (l,j)
    !	  fnetbi(j) = fnetbi(j) + directu(l,j) - direc(l,j)
    !	 END DO
  	END IF
      END DO
    END DO
 
   
    !  DO  i = 1, nir
    !	 tiru = tiru + firu(i)
    !  END DO
    !
    !	   Load fluxes into interface common block
    !
    !  DO j = 1, nlayer
    !	 fir_up_aerad(j) = fupbi(j)
    !	 fir_dn_aerad(j) = fdownbi(j)
    !  END DO
  
  END SUBROUTINE radtran
  
  SUBROUTINE oppr(ia,iz,ja,jz,m1,aot11)
    !
    !	  **************************************************************
    !	  *  Purpose		 :  CaLculates optical properties      *
    !	  *			    such as single scattering albedo,  *
    !	  *			    asymmetry parameter, etc.	       *
    !	  *			    This routine is case dependent and *
    !	  *			    wiLL have to be repLaced by the    *
    !	  *			    user.			       *
    !	  *  Subroutines Called  :  None			       *
    !	  *  Input		 :  PAH2O, RDH2O, CO2, O3, ETC         *
    !	  *  Output		 :  TAUL, W0, G0, OPD,         *
    !	  * ************************************************************
    !
    !INCLUDE 'globrad.h'
    !
    !	  W0(NWAVE,NLAYER) : SINGLE SCATTERING ALBEDO *** delta scaled ***
    !	  G0(NWAVE,NLAYER) : ASYMMETRY PARAMETER *** delta scaled ***
    !	  OPD(NWAVE,NLAYER): cumulative OPTICAL DEPTH *** delta scaled ***
    !	  SFL(NWAVE)	   : SOLAR FLUX
    !	 uW0(NWAVE,NLAYER)  : unscaled SINGLE SCATTERING ALBEDO
    !	 uG0(NWAVE,NLAYER)  : unscaled ASYMMETRY PARAMETER
    !	 uTAUL(NWAVE,NLAYER): unscaled OPTICAL DEPTH of layer
    !
    !	  ASSUME THAT P IS SAME ON ALL SIGMA LEVELS. IF PSURFACE
    !	  VARIES A LOT, THEN WE MAY NEED TO CALCULATE TAUO3,
    !	  TAUCO2, TAUO2 FOR EACH POINT.
    !
    !	  NOTE : THE TOP LAYER IS A DUMMY. IT CONTAINS A DEFINED GAS
    !		 AMOUNT. DIFFERENT MODELS WILL REQUIRE DIFFERENT
    !		 TREATMENT OF THIS.
    !	  CALCULATE TOTAL OPTICAL DEPTH INCLUDING GASES. THEN
    !	  GIVEN THE AEROSOL OPTICAL DEPTHS AND CLOUD OPTICAL DEPTHS,
    !	  CALCULATE FINAL OPTICAL PROPERTIES. WE USE A DELTA
    !	  TWO STREAM APPROACH TO FIND W0, SINGLE SCATTERING ALBEDO,
    !	  G0, ASYMMMETRY PARAMETER, TAUL, LAYER OPTICAL DEPTH,
    !	  OPD, CUMULATIVE OPTICAL DEPTH TO BASE OF LAYER.
    !
    USE mem_globrad, ONLY: nlayer,nwave,ngroup,nrad,xsecta, &
  			   rdqext,qscat,qbrqs,ntotal,nprob, &
  			   epsilon,g,ptop,p,q,nsolp,contnm, &
  			   uw0,ug0,ir,ngauss,gangle,ta,tb,  &
			   wa,wb,ga,gb,tia,tib,wia,wib,gia, &
			   gib,alpha,gama,caseE,caseW,  &          !kml2
          		   caseG,wave,imie               !kml2
    
    USE mem_aerad, ONLY: iprocopio
    USE carma_fastjx, ONLY: daer 
    use chem1_list, ONLY: PhotojMethod
    use mem_chem1,  only: CHEMISTRY

    IMPLICIT NONE
    
    INTEGER,INTENT(IN) :: ia,iz,ja,jz,m1
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: taua,taus,g01,wol
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: gol
    INTEGER :: i,i1,j1,kk,k
    INTEGER :: ig
    INTEGER :: iradgas
    INTEGER :: j
    INTEGER :: l
    REAL    :: cco((iz-ia+1)*(jz-ja+1))
    REAL    :: den 
    REAL    :: denom
    REAL    :: fo
    REAL    :: pcorr
    REAL    :: qcorr
    REAL    :: ttas
    REAL    :: tauh2o((iz-ia+1)*(jz-ja+1),ntotal,nlayer)
    REAL    :: utaul((iz-ia+1)*(jz-ja+1),ntotal,nlayer)
  !  REAL    :: uw0((iz-ia+1)*(jz-ja+1),ntotal,nlayer)
    REAL    :: wot((iz-ia+1)*(jz-ja+1),ntotal)
    REAL    :: got((iz-ia+1)*(jz-ja+1),ntotal)
    INTEGER,DIMENSION((iz-ia+1)*(jz-ja+1)):: idaot
    REAL :: aot11((iz-ia+1)*(jz-ja+1))         !kml2
    INTEGER :: ij,iend,jjj,in
    !lfr
    INTEGER,DIMENSION(4) :: nwl,nwu,npl,npu
    REAL,DIMENSION(4) :: wl
    REAL :: cm
    REAL,DIMENSION(4,(iz-ia+1)*(jz-ja+1),m1) :: daer_
    
!kmlnew
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: taucld,wcld,gcld 
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: wolc,woice,worain,gl,gice,grain
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: DENC
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: taucldlw,taucldice,taurain
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),nlayer) :: CORR,REFFI
    
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),nwave) :: rdqextnew,wonew,gonew
    real X_teste
  
     
    CORR=0.0
    DENC=0.0
    REFFI=0.0
    taucldlw=0.0
    taucldice=0.0
    taurain=0.0
    wolc=0.0
    woice=0.0
    worain=0.0
    gl=0.0
    gice=0.0
    grain=0.0
    taucld=0.0
    wcld=0.0
    gcld=0.0
    rdqextnew=0.0
    wonew=0.0
    gonew=0.0

!kmlnew    
       
    iend=(iz-ia+1)*(jz-ja+1)



      IF (iprocopio == 1 .and. imie == 1) THEN
!            aot11=0.1          !TMP KML2
 
   	    DO ij=1,iend
		idaot(ij) = MAX(MIN(INT(10*((ANINT(10.*aot11(ij))/10.)+0.1)/2.),9),1)

  	        DO  l = 1,nwave
	         rdqextnew(ij,l) = caseE(idaot(ij),l) 
		 wonew(ij,l)	 = caseW(idaot(ij),l) 
                 gonew(ij,l)	 = caseG(idaot(ij),l)
                 !if(l.eq.11) print*,'ext,wo,go=',rdqextnew(ij,l),wonew(ij,l),gonew(ij,l),ij
                END DO
            END DO

!TMP       
           taua=0.0
           DO j=1,nlayer
             DO ig = 1,ngroup
         	DO  i = 1,nrad
         	  DO  l = 1,ntotal
         	    DO ij=1,iend

         		taua(ij,l,j)=taua(ij,l,j)+rdqextnew(ij,nprob(l))*xsecta(i,ig)* &
         				   caer(ij,j,i,ig)
			tauaer(ij,l,j) = MAX(taua(ij,nprob(l),j),REAL(epsilon))		   
         		wol(ij,l,j)    = wonew(ij,nprob(l))
         		gol(ij,l,j)    = gonew(ij,nprob(l))
!               		if(tauaer(ij,l,j).gt.REAL(epsilon)) then
!			  print*,'ij,j,nprob=',ij,l,nprob(l)
!			  print*,'wo,go,ext,taua=',wol(ij,l,j),gol(ij,l,j),rdqextnew(ij,nprob(l))
!			  print*,'caer, taua=',taua(ij,l,j),caer(ij,j,i,ig)
!			endif
         	    END DO
         	  END DO
         	END DO
             END DO

           END DO
      ELSE

     	  taua=0.0
     	  taus=0.0
     	  g01=0.0
     	  DO j=1,nlayer
     	    DO ig = 1,ngroup
     	       DO  i = 1,nrad
     		 DO  l = 1,nwave
     		   DO ij=1,iend
     		     taua(ij,l,j)=taua(ij,l,j)+rdqext(i,ig,l)*xsecta(i,ig)* &
     				  caer(ij,j,i,ig)
     		     taus(ij,l,j)=taus(ij,l,j)+qscat(i,ig,l)*xsecta(i,ig)* &
     				  caer(ij,j,i,ig)
     		     g01(ij,l,j) =g01(ij,l,j) +qbrqs(i,ig,l)*xsecta(i,ig)* &
     				  caer(ij,j,i,ig)
     		   END DO
     		 END DO
     	       END DO
     	    END DO


     	    DO l= 1,ntotal
     	       DO ij=1,iend
     		 tauaer(ij,l,j) = MAX(taua(ij,nprob(l),j),REAL(epsilon))
!     		 tauaer(ij,l,j) =     taua(ij,nprob(l),j)
    		 wol(ij,l,j)	 = taus(ij,nprob(l),j)/tauaer(ij,l,j)
     		 ttas=1.0
     		 IF( wol(ij,l,j) /= 0. ) ttas = taus(ij,nprob(l),j)
     		 gol(ij,l,j)	= g01(ij,nprob(l),j)/ttas
     	       END DO
     	    END DO
     	  END DO
	  imie = 1
      END IF

!kmlnew 

!     imie = 1 

    DO j=1,nlayer
       DO ij=1,iend
       
       IF (xland_aerad(ij).ge..009) THEN
         REFFI(ij,j) =  7.0 * 1.e+3 * LWL_aerad(ij,j) + 5.5 
       ELSE 
         REFFI(ij,j) =  9.5 * 1.e+3 * LWL_aerad(ij,j) + 4.0 
       END IF
       
       CORR(ij,j) = 1.047 - 0.913e-4 * (tt(ij,j)-273.16) + 0.203e-3 * &
                   (tt(ij,j)-273.16) **2 - 0.106e-4 * (tt(ij,j)-273.16) **3
      
       CORR(ij,j) = MAX(CORR(ij,j),REAL((epsilon)))
              
	 
       END DO  
    END DO

   DO j=1,nlayer
      DO l= 1,ntotal
       DO ij=1,iend
        
	IF( j .eq. 1) taurain(ij,l,j) = 0.00018 * RAIN_aerad(ij) * 2000.0
	
	taucldlw(ij,l,j)= 1.e+3 * LWP_aerad(ij,j) *(ta(l)/REFFI(ij,j)+tb(l)/REFFI(ij,j)**2)

	worain(ij,l,j) = 1. - 0.45
	wolc(ij,l,j) = (1. - wa(l)) + wb(l) * REFFI(ij,j)
	gl(ij,l,j)   = ga(l) + gb(l) * REFFI(ij,j)
	grain(ij,l,j) = 0.95

	DENC(ij,l,j)      = 1. / (tia(l) + tib(l) * 1.0e+3 * IWL_aerad(ij,j))
	taucldice(ij,l,j) = CORR(ij,j) * 1.0e+3 * IWP_aerad(ij,j) * DENC(ij,l,j)
	
!srf-evitando bug - 0**0	
        if(wib(l) < 1.e-5) then
	   X_teste = 0.
	else
	   X_teste = (1.0e+3 * IWL_aerad(ij,j))**wib(l)
	endif

!	woice(ij,l,j) =( 1.0 -  wia(l) * (1.0e+3 * IWL_aerad(ij,j))**wib(l)) * & 
	woice(ij,l,j) =( 1.0 -  wia(l) *  X_teste                          ) * & 
	               ( 1.0 - gama(l) * (CORR(ij,j) - 1)/ CORR(ij,j) )
	
!srf-evitando bug - 0**0	
        if(gib(l) < 1.e-5) then
!	   X_teste = 0.
	   X_teste = 1.
	   
	else
	   X_teste = (1.0e+3 * IWL_aerad(ij,j))**gib(l)
	endif
!	gice(ij,l,j) = ( gia(l) * (1.0e+3 * IWL_aerad(ij,j))**gib(l) ) * &
	gice(ij,l,j) = ( gia(l) * X_teste                            ) * &
	               ( 1.0 + alpha(l) * (CORR(ij,j) - 1)/ CORR(ij,j) )
		       
	
	IF (l>=91 .AND. l<=113) THEN
	
	taucldlw(ij,l,j)= 1.0e+3 * LWP_aerad(ij,j) * ta(l) * exp(tb(l)* REFFI(ij,j))
	
	END IF
	
	IF (l>=114 .AND. l<=154) THEN
	
	taucldlw(ij,l,j)=  1.0e+3 * LWP_aerad(ij,j) * ( ta(l) + tb(l)* REFFI(ij,j))
	
	gl(ij,l,j) = 1. - ga(l) * exp( gb(l) * REFFI(ij,j))
	
	END IF
	
	taucld(ij,l,j)= taucldlw(ij,l,j) + taucldice(ij,l,j) + taurain(ij,l,j) !@@@@@@@@

!        if(LWL_aerad(ij,j).gt.0.) print*, 'Liquid', 1.0e+3 * LWL_aerad(ij,j),taucldlw(ij,l,j)
!	if(IWL_aerad(ij,j).gt.0.) print*, 'Ice', 1.0e+3 * IWL_aerad(ij,j),taucldice(ij,l,j)

        IF ( taucld(ij,l,j).gt.epsilon) THEN
	
	  wcld(ij,l,j) =  (wolc(ij,l,j) *  taucldlw(ij,l,j)  + &
	                  woice(ij,l,j) * taucldice(ij,l,j) + &
			  worain(ij,l,j) * taurain(ij,l,j)) &
	              / taucld(ij,l,j)
	  gcld(ij,l,j) = ( wolc(ij,l,j) *  taucldlw(ij,l,j)*   gl(ij,l,j) + &
	                woice(ij,l,j) * taucldice(ij,l,j)* gice(ij,l,j) + &
			worain(ij,l,j) * taurain(ij,l,j) * grain(ij,l,j)) &
		      / (wcld(ij,l,j) * taucld(ij,l,j))	     
        ELSE 
	  wcld(ij,l,j) = 1.0
	  gcld(ij,l,j) = 0.0
        ENDIF
       END DO                           
      END DO
    END DO
!kmlnew 
    
    iradgas = 1 !iradgas = 0: no gas in radiative xfer
    DO  j = 1,nlayer
      kk = MAX( 1, j-1 )
      !
      !   Bergstrom water vapor continuum fix:
      !
      !    <qcorr> is layer average water vapor mixing ratio
      !    <pcorr> is layer average pressure
      !
      !   For layer 0, calculate mixing ratio [g/g] from vapor column [g/cm^2]
      !   and average pressure [dyne/cm^2]
      !
      IF( j == 1 )THEN
        DO ij=1,iend
  	  qcorr = rdh2o(ij,1) * g / p_top(ij)
  	  pcorr = p_aerad(ij,1) / 2.
  	  cco(ij) = EXP(1800./t_aerad(ij,kk))*(qcorr*pcorr/2.87 + pcorr/4610.)
        END DO
      ELSE
        DO ij=1,iend
  	  qcorr = qv_aerad(ij,kk)
  	  pcorr = p_aerad(ij,kk)
  	  cco(ij) = EXP(1800./t_aerad(ij,kk))*(qcorr*pcorr/2.87 + pcorr/4610.)
        END DO
      END IF
  
      DO  l   = 1,ntotal
  	DO ij=1,iend
  	  IF (l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    tauh2o(ij,l,j) = rdh2o(ij,j)*pah2o(ij,l,j)
  	    !	  Bergstrom water vapor continuum fix (next two statements)
  	    IF( l > nsolp+30 .AND. l <= nsolp+36 ) THEN
  	      !kml	   if( L .GT. NSOLP+36 .AND. L .LE. NSOLP+42 ) then
  		  tauh2o(ij,l,j) = rdh2o(ij,j)*pah2o(ij,l,j)*cco(ij)
  	    ELSE
  		  tauh2o(ij,l,j) = tauh2o(ij,l,j)
  	    END IF
  	    IF (l > nsolp+36) tauh2o(ij,l,j) = tauh2o(ij,l,j) + &
  			    cco(ij)*rdh2o(ij,j)*contnm(l-nsolp)
     
  	    taul(ij,l,j)   = tauh2o(ij,l,j)+taugas(ij,l,j)+ &
  			      paray(ij,l,j)+tauaer(ij,l,j)+taucld(ij,l,j)
  
  	    IF (iradgas == 0) taul(ij,l,j) = tauaer(ij,l,j)
  	    IF( taul(ij,l,j) < epsilon ) taul(ij,l,j) = epsilon
  	  END IF
  	END DO
      END DO
  
      DO  l   = 1,ntotal
  	DO ij=1,iend

  	  IF (l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    utaul(ij,l,j)  = taul(ij,l,j)
  	    wot(ij,l)	   = (paray(ij,l,j)+tauaer(ij,l,j)*wol(ij,l,j)+  &
  			 taucld(ij,l,j)*wcld(ij,l,j))/taul(ij,l,j)
  	    IF (iradgas == 0) wot(ij,l) = wol(ij,l,j)
  	    wot(ij,l)	      = MIN(1.-REAL(epsilon),wot(ij,l))
  !	    uw0(ij,l,j)    = wot(ij)
  	    denom     = (paray(ij,l,j)+taucld(ij,l,j)*wcld(ij,l,j)+ &
  			 tauaer(ij,l,j)*wol(ij,l,j))
  	    !IF( denom <= epsilon ) denom = epsilon
  	    IF( denom > epsilon ) THEN
  	      got(ij,l) = ( wcld(ij,l,j)*gcld(ij,l,j)*taucld(ij,l,j) + &
  		      gol(ij,l,j)* wol(ij,l,j)*tauaer(ij,l,j) ) / denom
	      got(ij,l)=max(REAL(epsilon), got(ij,l))
  	    ELSE
  	      got(ij,l) = 0.
  	    END IF
  	    IF (iradgas == 0) got(ij,l) = gol(ij,l,j)
  	  END IF
  	END DO
      END DO

      DO  l   = 1,ntotal
  	DO ij=1,iend
  	  IF (l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    !ug0(l,j)	 = got(ij)
  	    fo        = got(ij,l)**2
  	    den       = 1.-wot(ij,l)*fo
  	    taul(ij,l,j)   = taul(ij,l,j) * den
  	    w0(ij,l,j)      = (1.-fo)*wot(ij,l)/den
  	    g0(ij,l,j)      = got(ij,l)/(1.+got(ij,l))
  	    opd(ij,l,j)    = 0.0
  	    opd(ij,l,j)    = opd(ij,l,kk)+taul(ij,l,j)
  	    uopd(ij,l,j)   = 0.0
  	    uopd(ij,l,j)   = uopd(ij,l,kk)+utaul(ij,l,j)
  	  END IF
  	END DO
      END DO
    END DO
  
    IF(ir_aerad == 1) THEN
      DO   j =   1,nlayer
  	DO  i =   1,ngauss
  	  DO  l =   1,ntotal
  	    DO ij=1,iend
  	      IF(l>=lls(ij) .AND. l<=lla(ij)) &
  		   y3(ij,l,i,j) = EXP(-taul(ij,l,j)/gangle(i))
  	    END DO
  	  END DO
  	END DO
      END DO
    END IF
    
    if (trim(PhotojMethod) == 'FAST-JX' .and. CHEMISTRY > 0)  then 
       return
       !LFR
       !Interpolating to obtain Opt. Depth aerosol
       !naer(1)=300nm - nwave=f(3)=296/f(4)=319    - nprob(3)/nprob(4) 
       !naer(2)=400nm - nwave=f(6)=365/f(7)=420    - nprob(6)/nprob(7)
       !naer(3)=600nm - nwave=f(13)=598/f(14)=660  - nprob(13)/nprob(14)
       !naer(4)=999nm - nwave=f(21)=926/f(22)=1005 - nprob(28)/nprob(32)
       nwl(1)=3;nwu(1)=4;wl(1)=.300
       nwl(2)=6;nwu(2)=7;wl(2)=.400
       nwl(3)=13;nwu(3)=14;wl(3)=.600
       nwl(4)=21;nwu(4)=22;wl(4)=.999
       npl(1)=3;npu(1)=4
       npl(2)=6;npu(2)=7
       npl(3)=13;npu(3)=14
       npl(4)=28;npu(4)=32
       DO k=1,m1
     	 DO ij=1,iend
     	   DO l=1,4
     	      cm=wl(l)/(wave(nwl(l))+wave(nwu(l)))
     	      daer_(l,ij,k)=(taucld(ij,npl(l),k)+taucld(ij,npu(l),k))*cm
     	   END DO
     	 END DO
       END DO
       !   
       CALL C_3d_4d(daer_,daer(:,:,ia:iz,ja:jz),m1,ia,iz,ja,jz,iend,4) 
       !LFR       
    endif 
  
  
  END SUBROUTINE oppr 
  
  SUBROUTINE oppr1(ia,iz,ja,jz,m1)
    !
    !	  **********************************************************
    !	  *  Purpose		 :  Calculate Planck Function and  *
    !	  *			    and its derivative at ground   *
    !	  *			    and at all altitudes.	   *
    !	  *  Subroutines Called  :  None			   *
    !	  *  Input		 :  TGRND, NLOW, WEIGHT 	   *
    !	  *  Output		 :  PTEMP, PTEMPG, SLOPE	   *
    !	  * ********************************************************
    !
    USE mem_globrad, ONLY: ntotal,tgrnd,nlow,nirp,plank,ltemp,nsolp, &
  			   weight,iblackbody_above,t_above, &
  			   nlayer,ncount
    
    IMPLICIT NONE
    
    INTEGER,INTENT(IN) :: ia,iz,ja,jz,m1
    INTEGER :: it1((iz-ia+1)*(jz-ja+1),nlayer)
    INTEGER :: itg((iz-ia+1)*(jz-ja+1))
    INTEGER :: itp((iz-ia+1)*(jz-ja+1))
    INTEGER :: j
    INTEGER :: kindex
    INTEGER :: l,i1,j1
    REAL :: pltemp1((iz-ia+1)*(jz-ja+1),ntotal)
    REAL :: ptemp2((iz-ia+1)*(jz-ja+1),ntotal,nlayer)
    INTEGER :: ij,iend,i
    
    iend=(iz-ia+1)*(jz-ja+1)
      !
    !	  **************************************
    !	  * CALCULATE PTEMP AND SLOPE	       *
    !	  **************************************
    !
    !	  CALCULATE THE WAVELENGTH DEPENDENT PLANK FUNCTION AT THE GROUND.
  
    DO ij=1,iend
      itg(ij)= ANINT(100.*t_surf(ij)) - nlow
!srf
!      if( itg(ij) < 0. .OR. itg(ij) > ncount) print*,'1-ITG=',itg(ij)
!srf      
    END DO
    DO i=1,nirp
      DO ij=1,iend
        pltemp1(ij,i)=plank(ltemp(i),itg(ij))
      END DO
    END DO
    DO  l =   nsolp+1,ntotal
      DO ij=1,iend
  	 ptempg(ij,l)=   pltemp1(ij,l-nsolp)*weight(l)
      END DO
    END DO
    !
    IF( iblackbody_above /= 0 )THEN
    !	    CALCULATE THE WAVELENGTH DEPENDENT PLANK FUNCTION AT THE TOP
    !	    OF THE MODEL.
      DO ij=1,iend
  	  itp(ij) = ANINT(100.*tabove_aerad(ij)) - nlow
!srf
!      if( itp(ij) < 0. .OR. itp(ij) > ncount) print*,'2-ITP=',itp(ij)
!srf      
      END DO
      DO i=1,nirp
        DO ij=1,iend  	 
	   pltemp1(ij,i)=plank(ltemp(i),itp(ij))
	   !CALL gather(nirp,pltemp1(ij,:),plank(1,itp(ij)),ltemp)
        END DO
      END DO
      DO  l =	nsolp+1,ntotal
  	DO ij=1,iend
  	    ptempt(ij,l)	=   pltemp1(ij,l-nsolp)*weight(l)
  	END DO
      END DO
  
    END IF
    !
    DO  j	     =   1,nlayer
      DO ij=1,iend
  	it1(ij,j) = ANINT(100.*tt(ij,j)) - nlow
!srf
!      if( it1(ij,j) < 0. .OR. it1(ij,j) > ncount) print*,'3-IT1=',it1(ij,j)
!srf      
      END DO
    END DO
    DO  j	     =   1,nlayer
      DO i=1,nirp
        DO ij=1,iend
           ptemp2(ij,i,j)=plank(ltemp(i),it1(ij,j))
        END DO
      END DO
    END DO
    
   ! kindex makes the top layer isothermal. using kindex, find
   ! plank function at bottom of each layer.
   ! note: if you force slope=0, then you have isothermal
   ! layers with tt(j) corresponding to average temperature
   ! of layer and tt(nlayer) should be set to tgrnd.
    DO  j	     =   1,nlayer
      kindex	      = MAX( 1, j-1 )
      DO  l	   = nsolp+1,ntotal
  	DO ij=1,iend
  	  ptemp(ij,l,j)   = ptemp2(ij,l-nsolp,j)*weight(l)
  	  slope(ij,l,j)   = (ptemp(ij,l,j)-ptemp(ij,l,kindex))/ &
  				 taul(ij,l,j)
  	  IF( taul(ij,l,j) <= 1.0E-6 ) slope(ij,l,j) = 0.
  	END DO
      END DO
    END DO
  
    !
  END SUBROUTINE oppr1
  
  SUBROUTINE twostr(m1,ia,iz,ja,jz)
    !
    !	 ******************************************************************
    !	 *  Purpose		:  Defines matrix properties and sets up  *
    !	 *			   matrix coefficients that do not depend *
    !	 *			   on zenith angle or temperature.	  *
    !	 *  Subroutines Called  :  None 				  *
    !	 *  Input		:  W0, G0				  *
    !	 *  Output		:  B1, B2, GAMI, ACON, EL1, AF, ETC	  *
    !	 * ****************************************************************
    !
  
    USE mem_globrad, ONLY: nsolp,sq3,tpi,nlayer,jn,jdble,irs,ntotal
  
    IMPLICIT NONE
  
    INTEGER,INTENT(IN) :: m1,ia,iz,ja,jz
    INTEGER	   :: j
    INTEGER	   :: jd
    INTEGER	   :: l
    REAL,PARAMETER :: two = 2.d0
    INTEGER :: ij,iend
    
    iend=(iz-ia+1)*(jz-ja+1)
      
    DO  l    =  1,ntotal !lls(i1,j1),lla(i1,j1)
      DO ij=1,iend  
  	IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	  IF( l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    IF(l <= nsolp ) THEN
  	      u1i(ij,l) = sq3
  	    ELSE
  	      u1i(ij,l) = two
  	    END IF
  	    !u1s(l)  =  tpi/u1i(l)
  	  END IF
  	END IF
      END DO
    END DO
    !
    !	   here we define layer properties following general scheme
    !	   of meador and weavor. then we set up layer properties
    !	   needed for matrix.
    !
    DO  j =  1,nlayer
      DO  l=  1,ntotal
  	DO  ij=  1,iend
  	  IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	    IF( l>=lls(ij) .AND. l<=lla(ij)) THEN
  	      !these are for two stream and hemispheric means
  	      b1(ij,l,j)   =  0.5*u1i(ij,l)*(2.-w0(ij,l,j)*(1. + g0(ij,l,j)))
  	      b2(ij,l,j)   =  0.5*u1i(ij,l)*w0(ij,l,j)*(1. - g0(ij,l,j))
  	      ak(ij,l,j)   =  SQRT(ABS(b1(ij,l,j)**2 - b2(ij,l,j)**2))
  	      gami(ij,l,j)  =  b2(ij,l,j)/(b1(ij,l,j) + ak(ij,l,j))
  	      ee1(ij,l,j)   =  EXP(-ak(ij,l,j)*taul(ij,l,j))
  	      el1(ij,l,j)   =  1.0 + gami(ij,l,j) *ee1(ij,l,j)
  	      em1(ij,l,j)   =  1.0 - gami(ij,l,j) * ee1(ij,l,j)
  	      el2(ij,l,j)   =  gami(ij,l,j) + ee1(ij,l,j)
  	      em2(ij,l,j)   =  gami(ij,l,j) - ee1(ij,l,j)
  	    END IF
  	  END IF
  	END DO
      END DO
    END DO
    !
    !	  we seek to solve ax(l-1)+bx(l)+ex(l+1) = d.
    !	  l=2n for even l, l=n+1 for odd l. the mean intensity (tmi/4pi)
    !	  and the net flux (fnet) are related to x's as noted in add.
    !	  first we set up the coefficients that are independent of solar
    !	  angle or temparature: a(i),b(i),e(i). d(i) is defined in add.
    !
    j=  0
    DO  jd=  2,jn,2
      j=  j + 1
      DO  l=  1,ntotal
  	DO  ij=  1,iend
  	  IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	    IF( l>=lls(ij) .AND. l<=lla(ij)) THEN
  	      !here are the even matrix elements
  	      af(ij,l,jd)   =  em1(ij,l,j+1)*el1(ij,l,j)- &
  				  em2(ij,l,j+1)*el2(ij,l,j)
  	      bf(ij,l,jd)   =  em1(ij,l,j+1)* em1(ij,l,j)- &
  				  em2(ij,l,j+1)*em2(ij,l,j)
  	      ef(ij,l,jd)   =  el1(ij,l,j+1)*em2(ij,l,j+1) - &
  				  el2(ij,l,j+1)*em1(ij,l,j+1)
  	      !here are the odd matrix elements except for the top.
  	      af(ij,l,jd+1) =  em1(ij,l,j)*el2(ij,l,j)- &
  				  el1(ij,l,j)*em2(ij,l,j)
  	      bf(ij,l,jd+1) =  el1(ij,l,j+1)*el1(ij,l,j) - &
  				  el2(ij,l,j+1)*el2(ij,l,j)
  	      ef(ij,l,jd+1) =  el2(ij,l,j)*em2(ij,l,j+1)- &
  				el1(ij,l,j)*em1(ij,l,j+1)
  	    END IF
  	  END IF
  	END DO
      END DO
    END DO
    !
    !	  HERE ARE THE TOP AND BOTTOM BOUNDARY CONDITIONS AS WELL AS THE
    !	  BEGINNING OF THE TRIDIAGONAL SOLUTION DEFINITIONS. I ASSUME
    !	  NO DIFFUSE RADIATION IS INCIDENT AT UPPER BOUNDARY.
    !
    DO  l=  1,ntotal
      DO  ij=  1,iend
  	IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	  IF( l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    af(ij,l,1)    = 0.0
  	    bf(ij,l,1)    = el1(ij,l,1)
  	    ef(ij,l,1)    = -em1(ij,l,1)
  	    af(ij,l,jdble) = el1(ij,l,nlayer)-rsfx(ij,l)*el2(ij,l,nlayer)
  	    bf(ij,l,jdble) = em1(ij,l,nlayer)-rsfx(ij,l)*em2(ij,l,nlayer)
  	    ef(ij,l,jdble) = 0.0
  	  END IF
  	END IF
      END DO
    END DO
  
  END SUBROUTINE twostr
  
  SUBROUTINE add(m1,ia,iz,ja,jz,cosz,m2,m3)
  
    USE mem_globrad, ONLY: isl,u0,nlayer,nsolp,sq3,sol,epsilon, &
  			   irs,ntotal,u1s,emis,pi,jn,tpi, &
  			   jdble,ndbl
    
    IMPLICIT NONE 
   
    !	  THIS SUBROUTINE FORMS THE MATRIX FOR THE MULTIPLE LAYERS AND
    !	  USES A TRIDIAGONAL ROUTINE TO FIND RADIATION IN THE ENTIRE
    !	  ATMOSPHERE.
   
    !	  ******************************
    !	  *   CALCULATIONS FOR SOLAR   *
    !	  ******************************
    INTEGER,INTENT(IN) :: ia,iz,ja,jz,m1,m2,m3
    REAL,INTENT(IN),DIMENSION((iz-ia+1)*(jz-ja+1)) :: cosz  
    INTEGER :: j,kk
    INTEGER :: jd
    INTEGER :: kindex
    INTEGER :: l
    REAL    :: b4
    REAL    :: c1
    REAL    :: c2
    REAL    :: cm1
    REAL    :: cp1
    REAL    :: du0
    REAL    :: x
    REAL    :: x2
    REAL    :: x3
    REAL    :: x4
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),nsolp,nlayer) :: direct,el3,ee3,cm
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),nsolp) :: sfcs
    REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,ndbl) :: df,as,ds,xk
    INTEGER :: ij,iend,i1,j1
    
    iend=(iz-ia+1)*(jz-ja+1)
  
    DO  j	 =  1,nlayer
      kk = MAX( 1, j-1 )
      DO  l    =  1,nsolp
  	DO ij=1,iend
  	  du0=1./cosz(ij)
  	  IF(isl_aerad(ij) /= 0)  THEN
  	    b3(ij,l,j)     =  0.5*(1.-sq3*g0(ij,l,j)*cosz(ij))
  	    b4         =  1. - b3(ij,l,j)
  	    x2         =  taul(ij,l,j)*du0
  	    ee3(ij,l,j)   =  EXP(-x2)
  	    x3         =  opd(ij,l,j)*du0
  	    el3(ij,l,j)   =  EXP(-x3)*sol(l)
  	    direct(ij,l,j) =  cosz(ij)*el3(ij,l,j)
  	    c1         =  b1(ij,l,j) - du0
  	    IF( ABS(c1) < epsilon ) c1 = SIGN(REAL(epsilon),c1)
  	    c2         =  ak(ij,l,j)*ak(ij,l,j) - du0*du0
  	    IF( ABS(c2) <= epsilon ) c2 = epsilon
  	    cp1        =  w0(ij,l,j)*(b3(ij,l,j)*c1+b4*b2(ij,l,j))/c2
  	    cpb(ij,l,j)    =  cp1 * el3(ij,l,j)
  	    IF( j /= 1 ) THEN
  	      x4 = el3(ij,l,kk)
  	    ELSE
  	      x4 = sol(l)
  	    END IF
  	    cp(ij,l,j)     =  cp1 * x4
  	    cm1        =  ( cp1*b2(ij,l,j) + w0(ij,l,j)*b4 )/c1
  	    cmb(ij,l,j)    =  cm1 * el3(ij,l,j)
  	    cm(ij,l,j)    =  cm1 * x4
  	  END IF
  	END DO
      END DO
    END DO
    !	     CALCULATE SFCS, THE SOURCE AT THE BOTTOM.
    DO  l=  1,nsolp
      DO ij=1,iend
  	IF(isl_aerad(ij) /= 0)  THEN
  	  sfcs(ij,l)=  direct(ij,l,nlayer) * rsfx(ij,l)
  	END IF
      END DO
    END DO
   
    !	  ******************************
    !	  * CALCULATIONS FOR INFRARED. *
    !	  ******************************
    DO  j= 1,nlayer
      DO  l = nsolp+1,ntotal
  	DO ij=1,iend
  	  IF(irs /= 0)  THEN
  	      kindex = MAX(1,j-1)
  	      b3(ij,l,j)     = 1.0/(b1(ij,l,j)+b2(ij,l,j))
  	      cp(ij,l,j)     = (ptemp(ij,l,kindex)+slope(ij,l,j)* &
  				   b3(ij,l,j))*(tpi/u1i(ij,l))
  	      cpb(ij,l,j)    = cp(ij,l,j) + slope(ij,l,j)* &
  				  taul(ij,l,j)*(tpi/u1i(ij,l))
  	      cm(ij,l,j)     = (ptemp(ij,l,kindex)-slope(ij,l,j)* &
  				   b3(ij,l,j))*(tpi/u1i(ij,l))
  	      cmb(ij,l,j)    = cm(ij,l,j) + slope(ij,l,j)* &
  				  taul(ij,l,j)*(tpi/u1i(ij,l))
  	      el3(ij,l,j)    = 0.0
  	      direct(ij,l,j) = 0.0
  	      ee3(ij,l,j)    = 0.0
  	  END IF
  	END DO
      END DO
    END DO
    
    DO  l= nsolp+1,ntotal
      DO ij=1,iend
  	IF(irs /= 0)  THEN
  	  sfcs(ij,l)= emis(l)*ptempg(ij,l)*pi
  	END IF
      END DO
    END DO
   
    j=  0
    DO  jd=  2,jn,2
     j=  j + 1
     DO  l=1,ntotal
       DO ij=1,iend
  	 IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	   IF(l>=lls(ij) .AND. l<=lla(ij)) THEN
  	 !	    HERE ARE THE EVEN MATRIX ELEMENTS
  	   df(ij,l,jd) = (cp(ij,l,j+1) - cpb(ij,l,j))*em1(ij,l,j+1) -  &
  		(cm(ij,l,j+1) - cmb(ij,l,j))*em2(ij,l,j+1)
  	 !	    HERE ARE THE ODD MATRIX ELEMENTS EXCEPT FOR THE TOP.
  	   df(ij,l,jd+1) =  el2(ij,l,j) * (cp(ij,l,j+1)-cpb(ij,l,j)) +  &
  		el1(ij,l,j) * (cmb(ij,l,j) - cm(ij,l,j+1))
  	    END IF
  	  END IF
  	END DO
      END DO
    END DO
   
    !	  HERE ARE THE TOP AND BOTTOM BOUNDARY CONDITIONS AS WELL AS THE
    !	  BEGINNING OF THE TRIDIAGONAL SOLUTION DEFINITIONS. I ASSUME NO
    !	  DIFFUSE RADIATION IS INCIDENT AT THE TOP.
    DO  l=1,ntotal
      DO ij=1,iend
  	IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	  IF(l>=lls(ij).AND. l<=lla(ij)) THEN
  	    df(ij,l,1)   = -cm(ij,l,1)
  	    df(ij,l,jdble) = sfcs(ij,l)+rsfx(ij,l)*cmb(ij,l,nlayer)- &
  			  cpb(ij,l,nlayer)
  	    ds(ij,l,jdble) = df(ij,l,jdble)/bf(ij,l,jdble)
  	    as(ij,l,jdble) = af(ij,l,jdble)/bf(ij,l,jdble)
  	  END IF
  	END IF
      END DO
    END DO
   
    !	  ********************************************
    !	  *	WE SOLVE THE TRIDIAGONAL EQUATIONS   *
    !	  ********************************************
   
    DO  j = 2, jdble
      DO  l=1,ntotal
  	DO ij=1,iend
  	  IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	    IF(l>=lls(ij) .AND. l<=lla(ij)) THEN
  	      x  = 1./(bf(ij,l,jdble+1-j) - ef(ij,l,jdble+1-j)* &
  			  as(ij,l,jdble+2-j))
  	      as(ij,l,jdble+1-j) = af(ij,l,jdble+1-j)*x
  	      ds(ij,l,jdble+1-j) = (df(ij,l,jdble+1-j) - &
  			  ef(ij,l,jdble+1-j) *ds(ij,l,jdble+2-j))*x
  	    END IF
  	  END IF
  	END DO
      END DO
    END DO
   
    DO  l=1,ntotal
      DO ij=1,iend
  	IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	  IF(l>=lls(ij) .AND. l<=lla(ij)) THEN
  	    xk(ij,l,1)    = ds(ij,l,1)
  	  END IF 
  	END IF
      END DO
    END DO
  
    DO  j	= 2, jdble
      DO  l=1,ntotal
  	DO ij=1,iend
  	  IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	    IF(l>=lls(ij) .AND. l<=lla(ij)) THEN
  	      xk(ij,l,j) = ds(ij,l,j) - as(ij,l,j)*xk(ij,l,j-1)
  	    END IF 
  	  END IF
  	END DO
      END DO
    END DO
   
    !  ***************************************************************
    !	  CALCULATE LAYER COEFFICIENTS, NET FLUX AND MEAN INTENSITY
    !  ***************************************************************
      
     DO j = 1,nlayer
       DO  l=1,ntotal
  	 DO ij=1,iend
  	   IF(isl_aerad(ij) /= 0 .OR. irs .NE. 0 ) THEN
  	     IF(l>=lls(ij) .AND. l<=lla(ij)) THEN
  	       ck1(ij,l,j)   = xk(ij,l,2*j-1)
  	       ck2(ij,l,j)   = xk(ij,l,2*j)
	       
  	       fnet(ij,l,j)  = ck1(ij,l,j)  *( el1(ij,l,j) &
  				 -el2(ij,l,j)) + ck2(ij,l,j) * &
  				    ( em1(ij,l,j)-em2(ij,l,j) ) + &
  				    cpb(ij,l,j) - cmb(ij,l,j) - direct(ij,l,j)
 	       tmi(ij,l,j)   =  el3(ij,l,j) + u1i(ij,l) *(ck1(ij,l,j)  *  &
  			  ( el1(ij,l,j) + el2(ij,l,j))   + ck2(ij,l,j) * &
  			  ( em1(ij,l,j)+em2(ij,l,j) ) +  cpb(ij,l,j) + &
  			  cmb(ij,l,j) )
  	     END IF 
  	   END IF
  	 END DO
       END DO
     END DO
   
    END SUBROUTINE add
    
    SUBROUTINE newflux1(m1,ia,iz,ja,jz)
      !
      !     **************************************************************
      !     *  Purpose  	   :  Calculate upward and downward	 *
      !     *			      intensities and fluxes using Gauss *
      !     *			      Quadrature angles and weights.	 *
      !     *  Subroutines Called  :  None				 *
      !     *  Input		   :  PTEMP, SLOPE, Y3, B3, EE1, EE2	 *
      !     *  Output		   :  DINTENT, UINTENT, DIREC, DIRECTU   *
      !     * ************************************************************
      !
      !INCLUDE 'globrad.h'
      USE mem_globrad, ONLY: ntotal,ngauss,nlayer,nsolp,tpi, &
  			     irs,gangle, &
  			     iblackbody_above, &
  			     gratio,gweight,emis
      
      IMPLICIT NONE
      
      INTEGER,INTENT(IN) :: m1,ia,iz,ja,jz
      INTEGER :: i
      INTEGER :: j
      INTEGER :: kindex
      INTEGER :: l
      INTEGER :: m
      REAL    :: ckm
      REAL    :: ckp
      REAL    :: x4
      REAL    :: ya
      REAL    :: yb
  	
      !
      REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,ngauss,nlayer) :: y1,y2,y4,y8
      REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,ngauss,nlayer) :: dintent,uintent
      REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: a1,a2,a3,a4,a7
      REAL,DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nlayer) :: y5
      INTEGER :: ij,iend,i1,j1
      
      iend=(iz-ia+1)*(jz-ja+1)
      !
    
      DO  j	  =  1,nlayer
  	kindex     = MAX( 1, j-1 )
  	DO   l  =  nsolp+1,ntotal
  	  DO ij=1,iend
  	    !HERE WE DO NO SCATTERING COEFFICIENTS
  	    a3(ij,l,j) =  ptemp(ij,l,kindex)*tpi
  	    a4(ij,l,j) =  tpi*slope(ij,l,j)
  	    a7(ij,l,j) =  a3(ij,l,j)
  	    y5(ij,l,j) =  a4(ij,l,j)*taul(ij,l,j)
  	  END DO
      END DO
      !HERE WE DO SCATTERING
      DO  l    =  nsolp+1,ntotal
  	DO ij=1,iend
  	  IF(irs /= 0) THEN
  	    x4         =  slope(ij,l,j)*(tpi*b3(ij,l,j)-(tpi/u1i(ij,l)))
  	    a1(ij,l,j) = u1i(ij,l) - ak(ij,l,j)
  	    a2(ij,l,j) = gami(ij,l,j)*(ak(ij,l,j)+u1i(ij,l))
  	    a3(ij,l,j) = a3(ij,l,j)+x4
  	    a7(ij,l,j) = a7(ij,l,j)-x4
  	  END IF
  	END DO
      END DO
    END DO
    !
    !	  CALCULATIONS FOR ALL GAUSS POINTS. HERE WE DO NO SCATTERING COEFFI
    !
    DO j=  1,nlayer
      DO i=  1,ngauss
  	DO  l =  nsolp+1,ntotal
  	  DO ij=1,iend
  	    y1(ij,l,i,j)  =  0.0
  	    y2(ij,l,i,j)  =  0.0
  	    y4(ij,l,i,j)  =  a7(ij,l,j) - a4(ij,l,j)*gangle(i)
  	    y8(ij,l,i,j)  =  a3(ij,l,j) + a4(ij,l,j)*gangle(i)
  	  END DO
  	END DO
  	!HERE WE DO SCATTERING
  	DO  l =  nsolp+1,ntotal
  	  DO ij=1,iend
  	    IF(irs /= 0) THEN
  	      ya=  a1(ij,l,j)*(y3(ij,l,i,j)-ee1(ij,l,j))/ &
  		      (ak(ij,l,j)*gangle(i)-1.)
  	      yb=  a2(ij,l,j)*(1.- ee1(ij,l,j)*y3(ij,l,i,j))/ &
  			(ak(ij,l,j)*gangle(i)+1.)
  	      ckp= ck1(ij,l,j)+ck2(ij,l,j)
  	      ckm= ck1(ij,l,j) -ck2(ij,l,j)
  	      y1(ij,l,i,j) =  ckp*yb+ckm*ya
  	      y2(ij,l,i,j) =  ckp*ya+ ckm*yb
  	    END IF
  	  END DO
  	END DO
      END DO
    END DO
    !
    DO  j	  =  1,nlayer
      DO   l	 =  nsolp+1,ntotal
  	DO ij=1,iend
  	  !tmid(ij,l,j) = 0.0
  	  !tmiu(ij,l,j) = 0.0
  	  direc(ij,l,j)     =  0.0
  	  directu(ij,l,j)   =  0.0
  	END DO
      END DO
    END DO
    !
    !	  DIREC IS DOWNWARD FLUX. DIRECTU IS UPWARD FLUX.
    !	  CALCULATE DINTENT THE DOWNWARD INTENSITY AND DIREC THE DOWNWARD FL
    !
    DO  i = 1,ngauss
      DO  l = nsolp+1,ntotal
  	DO ij=1,iend
  	  IF( iblackbody_above == 1 )THEN
  	    dintent(ij,l,i,1) = ptempt(ij,l)*y3(ij,l,i,1)*tpi + &
  		      y1(ij,l,i,1)+ (1.-y3(ij,l,i,1))*y4(ij,l,i,1)
  	  ELSE
  	    dintent(ij,l,i,1) = (1.-y3(ij,l,i,1))*y4(ij,l,i,1) + &
  		       y1(ij,l,i,1)
  	  END IF
  	  !tmid(ij,l,1) = tmid(ij,l,1)+dintent(ij,l,i,1)*gratio(i)
  	  direc(ij,l,1)= direc(ij,l,1)+dintent(ij,l,i,1)* gweight(i)
  	END DO
      END DO
    END DO
    !
    !	   DINTENT IS DOWNWARD INTENSITY * TPI. DIREC IS THE DOWNWARD FLUX.
    !
    DO j= 2,nlayer
      DO i = 1,ngauss
  	DO l = nsolp+1,ntotal
  	  DO ij=1,iend
  	    dintent(ij,l,i,j)  = dintent(ij,l,i,j-1)*y3(ij,l,i,j) + &
  			   y1(ij,l,i,j)+y5(ij,l,j)+  &
  			     (1.-y3(ij,l,i,j))*y4(ij,l,i,j)
  	    !tmid(ij,l,j)= tmid(ij,l,j)  +dintent(ij,l,i,j)*gratio(i)
  	    direc(ij,l,j)= direc(ij,l,j)+dintent(ij,l,i,j)*gweight(i)
  	  END DO
  	END DO
      END DO
    END DO
    !
    !	  UINTENT IS THE UPWARD INTENSITY * TPI. DIRECTU IS THE UPWARD FLUX.
    !	  ASSUME THAT THE REFLECTIVITY IS LAMBERT.
    !
    DO i =  1,ngauss
      DO l =  nsolp+1,ntotal
  	DO ij=1,iend
  	  uintent(ij,l,i,nlayer)  =  ptempg(ij,l)*emis(l) *tpi+2.* &
  			      rsfx(ij,l)*direc(ij,l,nlayer)
  	  !tmiu(ij,l,nlayer)=  tmiu(ij,l,nlayer)+ &
  	!		     uintent(ij,l,i,nlayer)*gratio(i)
  	  directu(ij,l,nlayer)    =  directu(ij,l,nlayer)+ &
  			      uintent(ij,l,i,nlayer)*gweight(i)
  	END DO
      END DO
    END DO
    !
    DO m= 2,nlayer
      j = nlayer-m+1
      DO i = 1,ngauss
  	DO l = nsolp+1,ntotal
  	  DO ij=1,iend
  	    uintent(ij,l,i,j)= (uintent(ij,l,i,j+1)-y5(ij,l,j+1)) * &
  			  y3(ij,l,i,j+1)+y2(ij,l,i,j+1)+ &
  			(1.-y3(ij,l,i,j+1))*y8(ij,l,i,j+1)
  	    !tmiu(ij,l,j) = tmiu(ij,l,j)+uintent(ij,l,i,j)*gratio(i)
  	    directu(ij,l,j) = directu(ij,l,j) + gweight(i)* &
  			       uintent(ij,l,i,j)
  	  END DO
  	END DO
      END DO
    END DO
    
  END SUBROUTINE newflux1
  
  
  SUBROUTINE plnk(e,t1,d)
    
    !	  ******************************************************
    !	  *  Purpose		 :  Calculate Planck Function  *
    !	  *  Subroutines Called  :  None		       *
    !	  *  Input		 :  WAVE, NCOUNT	       *
    !	  *  Output		 :  PLANK		       *
    !	  * ****************************************************
   
    !  THIS SUBROUTINE COMPUTES THE INTEGRAL OF THE PLANCK FUNCTION BETWEEN
    !  ZERO AND THE SPECIFIED VALUE OF LAMBDA.  THUS (USING XL AS LAMBDA)
    !  WE WANT TO INTEGRATE
    !  R = INTEGRAL(XL=0 TO XL=XLSPEC) ( C1*XL**-5* / (EXP(C2/XL*T)-1) )*DXL
    !  SUBSTITUTING U=C2/(XL*T), THE INTEGRAL BECOMES
    !  R = A CONSTANT TIMES INTEGRAL (USPEC TO INFINITY) OF
    !		 ( U**3 / (EXP(U) - 1) )*DU
    !  THE APPROXIMATIONS SHOWN HERE ARE ON PAGE 998 OF ABRAMOWITZ AND SEGUN
    !  UNDER THE HEADING OF DEBYE FUNCTIONS.  C2 IS THE PRODUCT OF PLANCK'S
    !  CONSTANT AND THE SPEED OF LIGHT DIVIDED BY BOLTZMANN'S CONSTANT.
    !  C2 = 14390 WHEN LAMBDA IS IN MICRONS.
    !  THE FACTOR 0.15399 IS THE RECIPROCAL OF SIX TIMES
    !  THE SUM OF (1/N**2) FOR ALL N FROM ONE TO INFINITY.  IT IS CHOSEN TO
    !  NORMALIZE THE INTEGRAL TO A MAXIMUM VALUE OF UNITY.
    !  RADIATION IN REAL UNITS IS OBTAINED BY MULTIPLYING THE INTEGRAL BY
    !  THE STEFAN-BOLTZMANN CONSTANT TIMES T**4.
    IMPLICIT NONE
   
    REAL				     :: e
    REAL, INTENT(IN)			     :: t1
    REAL, INTENT(OUT)			     :: d
    REAL :: am(5)
    REAL :: v1,a
    INTEGER :: m
   
    d		 =   0.0
    v1  	 =   e/t1
   
    IF (v1 <= 1.) THEN
      d 	=  1.0 - 0.15399*v1**3 *  &
  	  (1./3.-v1/8. + v1**2/60. - v1**4/5040. +  &
  	  v1**6/272160. - v1**8/13305600	 )
    END IF
   
    IF ( v1 > 1. .AND. v1 <= 50.) THEN
      DO  m   =  1,5
  	a	=  FLOAT(m)*v1
  	am(m)	=  0.15399 * EXP(-a)/m**4 * (((a+3.)*a+6.)*a+6.)
      END DO
   
      d 	 =  am(1)+am(2)+am(3)+am(4)+am(5)
    END IF
   
    d		  =  d*t1**4
   
  END SUBROUTINE plnk
  
  
  
  SUBROUTINE miess( ro, rfr, rfi, thetd, jx, qext, qscat,  &
  	  ctbrqs, r, re2, tmag2, wvno  )
   
    !
    ! **********************************************************************
    !	 THIS SUBROUTINE COMPUTES MIE SCATTERING BY A STRATIFIED SPHERE,
    !	 I.E. A PARTICLE CONSISTING OF A SPHERICAL CORE SURROUNDED BY A
    !	 SPHERICAL SHELL.  THE BASIC CODE USED WAS THAT DESCRIBED IN THE
    !	 REPORT: " SUBROUTINES FOR COMPUTING THE PARAMETERS OF THE
    !	 ELECTROMAGNETIC RADIATION SCATTERED BY A SPHERE " J.V. DAVE,
    !	 I B M SCIENTIFIC CENTER, PALO ALTO , CALIFORNIA.
    !	 REPORT NO. 320 - 3236 .. MAY 1968 .
    !
    !	 THE MODIFICATIONS FOR STRATIFIED SPHERES ARE DESCRIBED IN
    !	     TOON AND ACKERMAN, APPL. OPTICS, IN PRESS, 1981
    !
    !	 THE PARAMETERS IN THE CALLING STATEMENT ARE DEFINED AS FOLLOWS :
    !	   RO IS THE OUTER (SHELL) RADIUS;
    !	   R  IS THE CORE RADIUS;
    !	   RFR, RFI  ARE THE REAL AND IMAGINARY PARTS OF THE SHELL INDEX
    !	       OF REFRACTION IN THE FORM (RFR - I* RFI);
    !	   RE2, TMAG2  ARE THE INDEX PARTS FOR THE CORE;
    !	       ( WE ASSUME SPACE HAS UNIT INDEX. )
    !	   THETD(J): ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT
    !	       AND THE SCATTERED RADIATION.  THETD(J) IS< OR= 90.0
    !	       IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH
    !	       SUPPLEMENTARY VALUE, SEE COMMENTS BELOW ON ELTRMX;
    !	   JX: TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE
    !	       REQUIRED.  JX SHOULD NOT EXCEED IT UNLESS THE DIMENSIONS
    !	       STATEMENTS ARE APPROPRIATEDLY MODIFIED;
    !
    !	   THE DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN"LIGHT
    !	       SCATTERING BY SMALL PARTICLES,H.C.VAN DE HULST, JOHN WILEY '
    !	       SONS, INC., NEW YORK, 1957" .
    !	   QEXT: EFFICIENCY FACTOR FOR EXTINCTION,VAN DE HULST,P.14 ' 127.
    !	   QSCAT: EFFICIENCY FACTOR FOR SCATTERING,V.D. HULST,P.14 ' 127.
    !	   CTBRQS: AVERAGE(COSINE THETA) * QSCAT,VAN DE HULST,P.128
    !	   ELTRMX(I,J,K): ELEMENTS OF THE TRANSFORMATION MATRIX F,V.D.HULST
    !	       ,P.34,45 ' 125. I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1..
    !	       I = 3: ELEMENT S SUB 21.. I = 4: ELEMENT D SUB 21..
    !	   ELTRMX(I,J,1) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR
    !	       THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT
    !	       OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) ..
    !	   QBS IS THE BACK SCATTER CROSS SECTION.
    !
    !	   IT: IS THE DIMENSION OF THETD, ELTRMX, CSTHT, PI, TAU, SI2THT,
    !	       IT MUST CORRESPOND EXACTLY TO THE SECOND DIMENSION OF ELTRMX.
    !	   IACAP IS THE DIMENSION OF ACAP
    !	       IN THE ORIGINAL PROGRAM THE DIMENSION OF ACAP WAS 7000.
    !	       FOR CONSERVING SPACE THIS SHOULD BE NOT MUCH HIGHER THAN
    !	       THE VALUE, N=1.1*(NREAL**2 + NIMAG**2)**.5 * X + 1
    !	   WVNO: 2*PI / WAVELENGTH
    !
    !	 ALSO THE SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE O
    !	 DOWNWARD RECURRENCE RELATIONSHIP.
    !
    !	   TA(1): REAL PART OF WFN(1).  TA(2): IMAGINARY PART OF WFN(1).
    !	   TA(3): REAL PART OF WFN(2).  TA(4): IMAGINARY PART OF WFN(2).
    !	   TB(1): REAL PART OF FNA.	TB(2): IMAGINARY PART OF FNA.
    !	   TC(1): REAL PART OF FNB.	TC(2): IMAGINARY PART OF FNB.
    !	   TD(1): REAL PART OF FNAP.	TD(2): IMAGINARY PART OF FNAP.
    !	   TE(1): REAL PART OF FNBP.	TE(2): IMAGINARY PART OF FNBP.
    !	   FNAP, FNBP  ARE THE PRECEDING VALUES OF FNA, FNB RESPECTIVELY.
    ! **********************************************************************
    !
    !
    !	Include implicit declarations
    !
    !INCLUDE 'precision.h'
    !
    !
    !	Define dimensions of local arrays and arrays passed as arguments
    !
    IMPLICIT NONE
    
    INTEGER, PARAMETER :: iacap = 200000
    INTEGER, PARAMETER :: it = 1
    REAL, INTENT(IN)			     :: ro
    REAL, INTENT(IN)			     :: rfr
    REAL, INTENT(IN)			     :: rfi
    INTEGER, INTENT(IN) 		     :: jx
    REAL, INTENT(IN OUT)		     :: thetd(jx)
    REAL, INTENT(OUT)			     :: qext
    REAL, INTENT(OUT)			     :: qscat
    REAL, INTENT(OUT)			     :: ctbrqs
    REAL, INTENT(IN)			     :: r
    REAL, INTENT(IN)			     :: re2
    REAL, INTENT(IN)			     :: tmag2
    REAL, INTENT(IN)			     :: wvno
    !
    !
    !	Declare arguments passed as arrays
    !
   
    !
    !
    !	Declare local variables
    !
    DOUBLE PRECISION, PARAMETER :: epsilon_mie = 1.d-14
   
    DOUBLE COMPLEX fnap,   fnbp,   acap(iacap),  &
  	fna,	fnb,	rf,	  rrf, rrfx,   wm1,    fn1,	 fn2,  &
  	tc1,	tc2,	wfn(2),   z(4), k1,	k2,	k3,	  w(3,iacap),  &
  	rc,	u(8),	dh1, dh2,    dh4,    p24h24,   p24h21,  &
  	pstore, hstore, dummy,    dumsq
   
    DOUBLE PRECISION :: t(5), ta(4), tb(2), tc(2), td(2), te(2),  &
  	pi(3,it), tau(3,it), cstht(it), si2tht(it), eltrmx(4,it,2),  &
  	x, x1, x4, y1, y4, rx, sinx1, sinx4, cosx1, cosx4,  &
  	ey1, e2y1, ey4, ey1my4, ey1py4, aa, bb, cc, dd, denom,  &
  	realp, amagp, qbsr, qbsi, rmm
   
    INTEGER :: imag
    INTEGER :: iflag
    INTEGER :: nmx1
    INTEGER :: nmx2
    INTEGER :: n
    INTEGER :: nn
    INTEGER :: m
    INTEGER :: j
    INTEGER :: k
    INTEGER :: i
  
    EQUIVALENCE (fna,tb(1)),(fnb,tc(1)),(fnap,td(1)),(fnbp,te(1))
    !
    !
    !  Some compilers (e.g. absoft) don't the support imag(z) generic intrinsic function.
    !  In that situation, uncomment the following 2 lines to define stmt function to
    !  redefine imag() function to the alternative function.  I.e. change "alt_function"
    !  in following stmt function to the appropriate function that will return the
    !  imaginary part of a double complex for the compiler being used.
    !  For compilers that support imag(z), simply leave following 2 statments commented.
    !  The "c--alt_imag" comment prefix is designed to be recognized by the
    !  automated editing script create_dmiess, so please do not modify this prefix
    !  in the dmiess.f.template file.
    !  -bm  Nov-1999
    !
    DOUBLE COMPLEX z_dum_arg
   
    imag(z_dum_arg) = DIMAG(z_dum_arg)
    !
   
    !	    print*, RO, RFR, RFI, THETD, JX, QEXT, QSCAT, CTBRQS,
    !	  1		     R, RE2, TMAG2, WVNO
   
    !
    !	IF THE CORE IS SMALL SCATTERING IS COMPUTED FOR THE SHELL ONLY
    !
    iflag = 1
    IF ( r/ro < 1.d-6 )   iflag = 2
    IF ( jx <= it )   GO TO 20
    WRITE( *,7 )
    WRITE( *,6 )
    STOP 30
    20 rf =  CMPLX( rfr,  -rfi )
    rc =  CMPLX( re2, -tmag2 )
    x  =  ro * wvno
    k1 =  rc * wvno
    k2 =  rf * wvno
    k3 =  CMPLX( wvno, 0.0 )
    z(1) =  k2 * ro
    z(2) =  k3 * ro
    z(3) =  k1 * r
    z(4) =  k2 * r
    x1   =  REAL( z(1) )
    x4   =  REAL( z(4) )
    y1   =  DIMAG( z(1) )
    y4   =  DIMAG( z(4) )
    !	   print*,'Z(1)','Z(4)','x1','x4','y1','y4'
    !	   print*,Z(1),Z(4),x1,x4,y1,y4
   
    rrf  =  1.0 / rf
    rx   =  1.0 / x
    rrfx =  rrf * rx
    t(1) =  ( x**2 ) * ( rfr**2 + rfi**2 )
    t(1) =  SQRT( t(1) )
    nmx1 =  1.10 * t(1)
    !
    IF ( nmx1 <= iacap-1 )   GO TO 21
    WRITE(*,8)
    STOP 32
    21 nmx2 = t(1)
    IF ( nmx1 >  150 )   GO TO 22
    nmx1 = 150
    nmx2 = 135
    !
    22 acap( nmx1+1 )  =  ( 0.0,0.0 )
    IF ( iflag == 2 )	GO TO 26
    DO n = 1,3
      w( n,nmx1+1 )  =  ( 0.0,0.0 )
    END DO
    26 CONTINUE
    DO n = 1,nmx1
      nn = nmx1 - n + 1
      acap(nn) = (nn+1) * rrfx - 1.0 / ( (nn+1) * rrfx + acap(nn+1) )
      IF ( iflag == 2 )   GO TO 23
      DO m = 1,3
  	w( m,nn ) = (nn+1) / z(m+1)  - 1.0 / (  (nn+1) / z(m+1)  +  w( m,nn+1 )  )
      END DO
      23 CONTINUE
    END DO
    !
    DO    j = 1,jx
      IF ( thetd(j) < 0.0 )  thetd(j) =  ABS( thetd(j) )
      IF ( thetd(j) > 0.0 )  GO TO 24
      cstht(j)  = 1.0
      si2tht(j) = 0.0
      CYCLE
      24 IF ( thetd(j) >= 90.0 )  GO TO 25
      t(1)	=  ( 3.14159265359 * thetd(j) ) / 180.0
      cstht(j)  =  COS( t(1) )
      si2tht(j) =  1.0 - cstht(j)**2
      CYCLE
      25 IF ( thetd(j) > 90.0 )  GO TO 28
      cstht(j)  =  0.0
      si2tht(j) =  1.0
      CYCLE
      28 WRITE( *,5 )  thetd(j)
      WRITE( *,6 )
      STOP 34
    END DO
    !
    DO   j = 1,jx
      pi(1,j)  =  0.0
      pi(2,j)  =  1.0
      tau(1,j) =  0.0
      tau(2,j) =  cstht(j)
    END DO
    !
    ! INITIALIZATION OF HOMOGENEOUS SPHERE
    !
    t(1)   =  COS(x)
    t(2)   =  SIN(x)
    wm1    =  CMPLX( t(1),-t(2) )
    wfn(1) =  CMPLX( t(2), t(1) )
    ta(1)  =  t(2)
    ta(2)  =  t(1)
    wfn(2) =  rx * wfn(1) - wm1
    ta(3)  =  REAL(wfn(2))
    ta(4)  =  DIMAG(wfn(2))
    !	   print*,'WFN(2)','TA(3)','TA(4)'
    !	   print*,WFN(2),TA(3),TA(4)
   
    !
    IF ( iflag == 2 )	GO TO 560
    n = 1
    !
    ! INITIALIZATION PROCEDURE FOR STRATIFIED SPHERE BEGINS HERE
    !
    sinx1   =  SIN( x1 )
    sinx4   =  SIN( x4 )
    cosx1   =  COS( x1 )
    cosx4   =  COS( x4 )
    ey1     =  EXP( y1 )
    e2y1    =  ey1 * ey1
    ey4     =  EXP( y4 )
    ey1my4  =  EXP( y1 - y4 )
    ey1py4  =  ey1 * ey4
    ey1my4  =  EXP( y1 - y4 )
    aa  =  sinx4 * ( ey1py4 + ey1my4 )
    bb  =  cosx4 * ( ey1py4 - ey1my4 )
    cc  =  sinx1 * ( e2y1 + 1.0 )
    dd  =  cosx1 * ( e2y1 - 1.0 )
    denom   =  1.0  +  e2y1 * ( 4.0 * sinx1 * sinx1 - 2.0 + e2y1 )
    realp   =  ( aa * cc  +  bb * dd ) / denom
    amagp   =  ( bb * cc  -  aa * dd ) / denom
    dummy   =  CMPLX( realp, amagp )
    aa  =  sinx4 * sinx4 - 0.5
    bb  =  cosx4 * sinx4
    p24h24  =  0.5 + CMPLX( aa,bb ) * ey4 * ey4
    aa  =  sinx1 * sinx4  -  cosx1 * cosx4
    bb  =  sinx1 * cosx4  +  cosx1 * sinx4
    cc  =  sinx1 * sinx4  +  cosx1 * cosx4
    dd  = -sinx1 * cosx4  +  cosx1 * sinx4
    p24h21  =  0.5 * CMPLX( aa,bb ) * ey1 * ey4  + 0.5 * CMPLX( cc,dd ) * ey1my4
    dh4  =  z(4) / ( 1.0 + ( 0.0,1.0 ) * z(4) )  -  1.0 / z(4)
    dh1  =  z(1) / ( 1.0 + ( 0.0,1.0 ) * z(1) )  -  1.0 / z(1)
    dh2  =  z(2) / ( 1.0 + ( 0.0,1.0 ) * z(2) )  -  1.0 / z(2)
    pstore  =  ( dh4 + n / z(4) )  *  ( w(3,n) + n / z(4) )
    p24h24  =  p24h24 / pstore
    hstore  =  ( dh1 + n / z(1) )  *  ( w(3,n) + n / z(4) )
    p24h21  =  p24h21 / hstore
    pstore  =  ( acap(n) + n / z(1) )  /  ( w(3,n) + n / z(4) )
    dummy   =  dummy * pstore
    dumsq   =  dummy * dummy
    !
    ! NOTE:  THE DEFINITIONS OF U(I) IN THIS PROGRAM ARE NOT THE SAME AS
    !	     THE USUBI DEFINED IN THE ARTICLE BY TOON AND ACKERMAN.  THE
    !	     CORRESPONDING TERMS ARE:
    !	       USUB1 = U(1)			  USUB2 = U(5)
    !	       USUB3 = U(7)			  USUB4 = DUMSQ
    !	       USUB5 = U(2)			  USUB6 = U(3)
    !	       USUB7 = U(6)			  USUB8 = U(4)
    !	       RATIO OF SPHERICAL BESSEL FTN TO SPHERICAL HENKAL FTN = U(8)
    !
    u(1) =  k3 * acap(n)  -  k2 * w(1,n)
    u(2) =  k3 * acap(n)  -  k2 * dh2
    u(3) =  k2 * acap(n)  -  k3 * w(1,n)
    u(4) =  k2 * acap(n)  -  k3 * dh2
    u(5) =  k1 *  w(3,n)  -  k2 * w(2,n)
    u(6) =  k2 *  w(3,n)  -  k1 * w(2,n)
    u(7) =  ( 0.0,-1.0 )  *  ( dummy * p24h21 - p24h24 )
    u(8) =  ta(3) / wfn(2)
    !
    fna  =  u(8) * ( u(1)*u(5)*u(7)  +  k1*u(1)  -  dumsq*k3*u(5) ) /  &
  	( u(2)*u(5)*u(7)  +  k1*u(2)  -  dumsq*k3*u(5) )
    fnb  =  u(8) * ( u(3)*u(6)*u(7)  +  k2*u(3)  -  dumsq*k2*u(6) ) /  &
  	( u(4)*u(6)*u(7)  +  k2*u(4)  -  dumsq*k2*u(6) )
    GO TO 561
    560 tc1  =  acap(1) * rrf  +  rx
    tc2  =  acap(1) * rf   +  rx
    fna  =  ( tc1 * ta(3)  -  ta(1) ) / ( tc1 * wfn(2)  -  wfn(1) )
    fnb  =  ( tc2 * ta(3)  -  ta(1) ) / ( tc2 * wfn(2)  -  wfn(1) )
    !
    561 CONTINUE
    fnap = fna
    fnbp = fnb
    t(1) = 1.50
    !
    !	 FROM HERE TO THE STATMENT NUMBER 90, ELTRMX(I,J,K) HAS
    !	 FOLLOWING MEANING:
    !	 ELTRMX(1,J,K): REAL PART OF THE FIRST COMPLEX AMPLITUDE.
    !	 ELTRMX(2,J,K): IMAGINARY PART OF THE FIRST COMPLEX AMPLITUDE.
    !	 ELTRMX(3,J,K): REAL PART OF THE SECOND COMPLEX AMPLITUDE.
    !	 ELTRMX(4,J,K): IMAGINARY PART OF THE SECOND COMPLEX AMPLITUDE.
    !	 K = 1 : FOR THETD(J) AND K = 2 : FOR 180.0 - THETD(J)
    !	 DEFINITION OF THE COMPLEX AMPLITUDE: VAN DE HULST,P.125.
    !
    tb(1) = t(1) * tb(1)
    tb(2) = t(1) * tb(2)
    tc(1) = t(1) * tc(1)
    tc(2) = t(1) * tc(2)
    DO  j = 1,jx
      eltrmx(1,j,1) = tb(1) * pi(2,j) + tc(1) * tau(2,j)
      eltrmx(2,j,1) = tb(2) * pi(2,j) + tc(2) * tau(2,j)
      eltrmx(3,j,1) = tc(1) * pi(2,j) + tb(1) * tau(2,j)
      eltrmx(4,j,1) = tc(2) * pi(2,j) + tb(2) * tau(2,j)
      eltrmx(1,j,2) = tb(1) * pi(2,j) - tc(1) * tau(2,j)
      eltrmx(2,j,2) = tb(2) * pi(2,j) - tc(2) * tau(2,j)
      eltrmx(3,j,2) = tc(1) * pi(2,j) - tb(1) * tau(2,j)
      eltrmx(4,j,2) = tc(2) * pi(2,j) - tb(2) * tau(2,j)
    END DO
    !
    qext   = 2.0 * ( tb(1) + tc(1))
    qscat  = ( tb(1)**2 + tb(2)**2 + tc(1)**2 + tc(2)**2 ) / 0.75
    ctbrqs = 0.0
    qbsr   = -2.0*(tc(1) - tb(1))
    qbsi   = -2.0*(tc(2) - tb(2))
    rmm    = -1.0
    n = 2
    65 t(1) = 2*n - 1
    t(2) =   n - 1
    t(3) = 2*n + 1
    DO   j = 1,jx
      pi(3,j)  = ( t(1) * pi(2,j) * cstht(j) - n * pi(1,j) ) / t(2)
      tau(3,j) = cstht(j) * ( pi(3,j) - pi(1,j) )  -  &
  	  t(1) * si2tht(j) * pi(2,j)  +  tau(1,j)
    END DO
    !
    ! HERE SET UP HOMOGENEOUS SPHERE
    !
    wm1    =  wfn(1)
    wfn(1) =  wfn(2)
    ta(1)  =  REAL(wfn(1))
    ta(2)  =  DIMAG(wfn(1))
    ta(4)  =  DIMAG(wfn(2))
    wfn(2) =  t(1) * rx * wfn(1)  -  wm1
    ta(3)  =  REAL(wfn(2))
   
    !	   print*,'WFN(1)','TA(1)','TA(2)'
    !	   print*,WFN(1),TA(1),TA(2)
    !
    IF ( iflag == 2 )	GO TO 1000
    !
    ! HERE SET UP STRATIFIED SPHERE
    !
    dh2  =  - n / z(2)  +  1.0 / ( n / z(2) - dh2 )
    dh4  =  - n / z(4)  +  1.0 / ( n / z(4) - dh4 )
    dh1  =  - n / z(1)  +  1.0 / ( n / z(1) - dh1 )
    pstore  =  ( dh4 + n / z(4) )  *  ( w(3,n) + n / z(4) )
    p24h24  =  p24h24 / pstore
    hstore  =  ( dh1 + n / z(1) )  *  ( w(3,n) + n / z(4) )
    p24h21  =  p24h21 / hstore
    pstore  =  ( acap(n) + n / z(1) )  /  ( w(3,n) + n / z(4) )
    dummy   =  dummy * pstore
    dumsq   =  dummy * dummy
    !
    u(1) =  k3 * acap(n)  -  k2 * w(1,n)
    u(2) =  k3 * acap(n)  -  k2 * dh2
    u(3) =  k2 * acap(n)  -  k3 * w(1,n)
    u(4) =  k2 * acap(n)  -  k3 * dh2
    u(5) =  k1 *  w(3,n)  -  k2 * w(2,n)
    u(6) =  k2 *  w(3,n)  -  k1 * w(2,n)
    u(7) =  ( 0.0,-1.0 )  *  ( dummy * p24h21 - p24h24 )
    u(8) =  ta(3) / wfn(2)
    !
    fna  =  u(8) * ( u(1)*u(5)*u(7)  +  k1*u(1)  -  dumsq*k3*u(5) ) /  &
  	( u(2)*u(5)*u(7)  +  k1*u(2)  -  dumsq*k3*u(5) )
    fnb  =  u(8) * ( u(3)*u(6)*u(7)  +  k2*u(3)  -  dumsq*k2*u(6) ) /  &
  	( u(4)*u(6)*u(7)  +  k2*u(4)  -  dumsq*k2*u(6) )
    !
    1000 CONTINUE
    tc1  =  acap(n) * rrf  +  n * rx
    tc2  =  acap(n) * rf   +  n * rx
    fn1  =  ( tc1 * ta(3)  -  ta(1) ) /  ( tc1 * wfn(2) - wfn(1) )
    fn2  =  ( tc2 * ta(3)  -  ta(1) ) /  ( tc2 * wfn(2) - wfn(1) )
    m	 =  wvno * r
    IF ( n < m )   GO TO 1002
    IF ( iflag == 2 )	GO TO 1001
    IF (ABS((fn1-fna)/fn1) < epsilon_mie .AND.  &
  	ABS(  ( fn2-fnb ) / fn2  ) < epsilon_mie  ) iflag = 2
    IF ( iflag == 1 )	GO TO 1002
    1001 fna  =  fn1
    fnb  =  fn2
    !
    1002 CONTINUE
    t(5)  =  n
    t(4)  =  t(1) / ( t(5) * t(2) )
    t(2)  =  (  t(2) * ( t(5) + 1.0 )  ) / t(5)
    !
    ctbrqs  =  ctbrqs  +  t(2) * ( td(1) * tb(1)  +  td(2) * tb(2)  &
  	+	    te(1) * tc(1)  +  te(2) * tc(2) )  &
  	+  t(4) * ( td(1) * te(1)  +  td(2) * te(2) )
    qext    =	qext  +  t(3) * ( tb(1) + tc(1) )
    !	  $	   T(3), TB(1), TC(1), QEXT
    t(4)    =  tb(1)**2 + tb(2)**2 + tc(1)**2 + tc(2)**2
    qscat   =  qscat  +  t(3) * t(4)
    rmm     =  -rmm
    qbsr    =  qbsr + t(3)*rmm*(tc(1) - tb(1))
    qbsi    =  qbsi + t(3)*rmm*(tc(2) - tb(2))
    !
    t(2)    =  n * (n+1)
    t(1)    =  t(3) / t(2)
    k = (n/2)*2
    DO  j = 1,jx
      eltrmx(1,j,1) = eltrmx(1,j,1)+t(1)*(tb(1)*pi(3,j)+tc(1)*tau(3,j))
      eltrmx(2,j,1) = eltrmx(2,j,1)+t(1)*(tb(2)*pi(3,j)+tc(2)*tau(3,j))
      eltrmx(3,j,1) = eltrmx(3,j,1)+t(1)*(tc(1)*pi(3,j)+tb(1)*tau(3,j))
      eltrmx(4,j,1) = eltrmx(4,j,1)+t(1)*(tc(2)*pi(3,j)+tb(2)*tau(3,j))
      IF ( k == n )  THEN
  	eltrmx(1,j,2) =eltrmx(1,j,2)+t(1)*(-tb(1)*pi(3,j)+tc(1)*tau(3,j))
  	eltrmx(2,j,2) =eltrmx(2,j,2)+t(1)*(-tb(2)*pi(3,j)+tc(2)*tau(3,j))
  	eltrmx(3,j,2) =eltrmx(3,j,2)+t(1)*(-tc(1)*pi(3,j)+tb(1)*tau(3,j))
  	eltrmx(4,j,2) =eltrmx(4,j,2)+t(1)*(-tc(2)*pi(3,j)+tb(2)*tau(3,j))
      ELSE
  	eltrmx(1,j,2) = eltrmx(1,j,2)+t(1)*(tb(1)*pi(3,j)-tc(1)*tau(3,j))
  	eltrmx(2,j,2) = eltrmx(2,j,2)+t(1)*(tb(2)*pi(3,j)-tc(2)*tau(3,j))
  	eltrmx(3,j,2) = eltrmx(3,j,2)+t(1)*(tc(1)*pi(3,j)-tb(1)*tau(3,j))
  	eltrmx(4,j,2) = eltrmx(4,j,2)+t(1)*(tc(2)*pi(3,j)-tb(2)*tau(3,j))
      END IF
    END DO
    !
    IF ( t(4) < epsilon_mie )	GO TO 100
    n = n + 1
    DO  j = 1,jx
      pi(1,j)	=   pi(2,j)
      pi(2,j)	=   pi(3,j)
      tau(1,j)  =  tau(2,j)
      tau(2,j)  =  tau(3,j)
    END DO
    fnap  =  fna
    fnbp  =  fnb
    IF ( n <= nmx2 )   GO TO 65
    !	      print*,N,NMX2
    WRITE( *,8 )
    STOP 36
    100 CONTINUE
    DO j = 1,jx
      DO k = 1,2
  	DO    i= 1,4
  	  t(i)  =  eltrmx(i,j,k)
  	END DO
  	eltrmx(2,j,k)  =      t(1)**2  +  t(2)**2
  	eltrmx(1,j,k)  =      t(3)**2  +  t(4)**2
  	eltrmx(3,j,k)  =  t(1) * t(3)  +  t(2) * t(4)
  	eltrmx(4,j,k)  =  t(2) * t(3)  -  t(4) * t(1)
      END DO
    END DO
    t(1)    =	 2.0 * rx**2
    qext    =	qext * t(1)
    qscat   =  qscat * t(1)
    ctbrqs  =  2.0 * ctbrqs * t(1)
   
   
   
    !	   RO IS THE OUTER (SHELL) RADIUS;
    !	   R  IS THE CORE RADIUS;
    !	   RFR, RFI  ARE THE REAL AND IMAGINARY PARTS OF THE SHELL INDEX
    !	       OF REFRACTION IN THE FORM (RFR - I* RFI);
    !	   RE2, TMAG2  ARE THE INDEX PARTS FOR THE CORE;
   
    !
    ! QBS IS THE BACK SCATTER CROSS SECTION
    !
    !	   PIG   = ACOS(-1.0)
    !	   RXP4  = RX*RX/(4.0*PIG)
    !	   QBS   = RXP4*(QBSR**2 + QBSI**2)
    !
    5  FORMAT( 10X,' THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN 90.0 DEGREES. IT IS ', e15.4 )
    6  FORMAT( // 10X, 'PLEASE READ COMMENTS.' // )
    7  FORMAT( // 10X, 'THE VALUE OF THE ARGUMENT JX IS GREATER THAN IT'//)
    8  FORMAT( // 10X, 'THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST GET DETAILED OUTPUT AND MODIFY SUBROUTINE' // )
    !
  END SUBROUTINE miess
  
  
  SUBROUTINE  radtran_to_rams(m1,m2,m3,fthrl,rlong,fthrs,rshort,aotr,ia,iz,ja,jz,mynum)
  
    USE mem_grid   , ONLY: nzpmax	  !INTENT(IN)
    USE mem_globrad, ONLY: nwave,ntotal,nprob
    !USE carma_fastjx, ONLY: odcld   
    
    IMPLICIT NONE
  
    INTEGER,INTENT(IN)  	       :: m1,m2,m3,ia,iz,ja,jz,mynum
    REAL,INTENT(OUT)  ,DIMENSION((iz-ia+1)*(jz-ja+1))           :: rshort
    REAL,INTENT(OUT)  ,DIMENSION((iz-ia+1)*(jz-ja+1))           :: rlong
    REAL,INTENT(OUT)  ,DIMENSION((iz-ia+1)*(jz-ja+1),nwave)     :: aotr
    REAL,INTENT(INOUT),DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax)	:: fthrl
    REAL,INTENT(INOUT),DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax)	:: fthrs
    INTEGER :: ij,iend
  
    !Local
  
    REAL, DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax)         :: duml
    REAL, DIMENSION((iz-ia+1)*(jz-ja+1),nzpmax)         :: dums
    REAL, DIMENSION((iz-ia+1)*(jz-ja+1),nwave)          :: dumaot
    REAL, DIMENSION((iz-ia+1)*(jz-ja+1),ntotal,nzpmax)  :: dum2aot
    
  
    INTEGER :: k,j1,i1
    INTEGER :: k1
    INTEGER :: kr
    INTEGER :: l
    INTEGER :: nzz
  
    iend=(iz-ia+1)*(jz-ja+1)
    !nzz = Vertical level number
    nzz = m1 - 1
    aotr=0.0
       
!srf- otimizado
!    DO k = 1,nzz
!       !  Reverse the vertical index when in cartesian coordinates
!      !
!       k1  = nzz + 1 - k
!       DO ij=1,iend
!	duml(ij,k1) = heati_aerad(ij,k)
!	dums(ij,k1) = heats_aerad(ij,k)
!       END DO
!     END DO
!  
!    ! Transfer values from CARMA grid to BRAMS grid
!     DO k=1,m1-1
!       kr = k + 1     ! K level in CARMA grid corresponds to K+1 level in BRAMS grid
!       DO ij=1,iend
!	 fthrl(ij,kr) = duml(ij,k)
!	 fthrs(ij,kr) = dums(ij,k)
!       END DO
!     END DO
   
    ! reverse the vertical and transfer values from CARMA grid to BRAMS grid
    DO k=2,m1
       kr = nzz+2- k 
       DO ij=1,iend
  	 fthrl(ij,k) = heati_aerad(ij,kr)
  	 fthrs(ij,k) = heats_aerad(ij,kr)
!	print*,k,fthrl(ij,k),fthrs(ij,k)
       END DO
    END DO
!srf


     DO ij=1,iend
  	rshort(ij) = solnet(ij)  ! total absorvido pela superficie
        rlong(ij)  = xirdown(ij) ! downward longwave  na superficie
        !if (rlong(ij).lt.10.) print*,'RADTRAN_TO_RAMS!!!',ij,rlong(ij)
    END DO
    
    dumaot=0.0
    !  DATA WAVE / 0.256, 0.280, 0.296, 0.319, 0.335, 0.365, 0.420,
    ! 2   0.690, 0.762, 0.719, 0.813, 0.862, 0.926, 1.005, 1.111,
    ! 3   1.333, 1.562, 1.770, 2.051, 2.210, 2.584, 3.284, 3.809,
    ! 4   4.292, 4.546, 4.878, 5.128, 5.405, 5.714, 6.061, 6.452,
    ! 5   6.897, 7.407, 8.333, 9.009, 10.309, 12.500, 13.889,	 
    ! 6   16.667,20.000, 26.316, 35.714, 62.50  		/

    DO l=1,ntotal
       DO k=1,m1				  
  	 DO ij=1,iend
           dum2aot(ij,nprob(l),k)=tauaer(ij,nprob(l),k)
  	 END DO 				  
       END DO					  
    END DO

    DO l=1,nwave
       DO k=1,m1				  
  	 DO ij=1,iend
  	   dumaot(ij,l) = dumaot(ij,l) + dum2aot(ij,l,k)   
!  	   dumaot(ij,l) = dumaot(ij,l) + tauaer(ij,l,k)   
  	 END DO 				  
       END DO					  
    END DO
    
    DO l=1,nwave
      DO ij=1,iend
  	aotr(ij,l)= dumaot(ij,l)
      END DO
    END DO

  END SUBROUTINE radtran_to_rams

!kmlnew  
  SUBROUTINE radcomp_carma(m1,m2,m3,ia,iz,ja,jz,solfac  &
       ,theta,pi0,pp,rv,RAIN,LWL,IWL,dn0,rtp,fthrd  &
       ,rtgt,f13t,f23t,glat,glon,rshort,rlong,albedt,cosz,rlongup  &
       ,mynum,fmapt,patch_area,npat)
!kmlnew      
       USE mem_carma, ONLY: carma
       USE mem_grid , ONLY: ngrid

       ! For specific optimization depending the type of machine
       use machine_arq, only: machine ! INTENT(IN)

       INTEGER,INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz,mynum,npat
  
       REAL,INTENT(IN)    :: solfac
       REAL,INTENT(IN)    :: theta(m1,m2,m3)
       REAL,INTENT(IN)    :: pi0(m1,m2,m3)
       REAL,INTENT(IN)    :: pp(m1,m2,m3)
       REAL,INTENT(IN)    :: rv(m1,m2,m3)
!kmlnew
       REAL,INTENT(IN)    :: LWL(m1,m2,m3)
       REAL,INTENT(IN)    :: IWL(m1,m2,m3)
       REAL,INTENT(IN)    :: RAIN(m2,m3)
       REAL,INTENT(IN)    :: patch_area(m2,m3,npat)
!kmlnew      
       REAL,INTENT(IN)    :: dn0(m1,m2,m3)
       REAL,INTENT(IN)    :: rtp(m1,m2,m3)
       REAL,INTENT(IN)    :: rtgt(m2,m3)
       REAL,INTENT(IN)    :: f13t(m2,m3)
       REAL,INTENT(IN)    :: f23t(m2,m3)
       REAL,INTENT(IN)    :: glat(m2,m3)
       REAL,INTENT(IN)    :: glon(m2,m3)
       REAL,INTENT(IN)    :: cosz(m2,m3) 
       REAL,INTENT(IN)    :: albedt(m2,m3)
       REAL,INTENT(IN)    :: fmapt(m2,m3)
!       REAL,INTENT(IN)    :: pm(m1,m2,m3) 
  
       REAL,INTENT(INOUT) :: rshort(m2,m3)
       REAL,INTENT(INOUT) :: rlong(m2,m3)
       REAL,INTENT(INOUT) :: fthrd(m1,m2,m3)
       
       REAL,INTENT(INOUT) :: rlongup(m2,m3)
       !kmlnew
       REAL :: xland(m2,m3)
       
       INTEGER :: ia1,iz1,ja1,jz1,ii,jj,iend,ipat
       
       iend=(iz-ia+1)*(jz-ja+1)

       ! For specific optimization
       if (machine==1) then
          !-sx6 
          p_isize=16;p_jsize=16
       elseif(machine==0) then
          !-cluster ! Generic IA32
          p_isize=2;p_jsize=2
       endif
       
!srf      DO ipat= 1,2
       DO ii=ia,iz,p_isize
  	 DO jj=ja,jz,p_jsize
  	   ia1=min(ii,iz);iz1=min(ii+p_isize-1,iz)
  	   ja1=min(jj,jz);jz1=min(jj+p_jsize-1,jz)
  	 END DO
       END DO
!srf      END DO           

       DO ii=ia,iz
  	 DO jj=ja,jz
	   xland(ii,jj) = patch_area(ii,jj,2)
  	 END DO
       END DO

      
       DO ii=ia,iz,p_isize
  	 DO jj=ja,jz,p_jsize
  	   ia1=min(ii,iz);iz1=min(ii+p_isize-1,iz)
  	   ja1=min(jj,jz);jz1=min(jj+p_jsize-1,jz)

  	   CALL radcarma(m1,m2,m3,ia1,iz1,ja1,jz1,solfac  &
  			    ,theta,pi0,pp,rv,RAIN,LWL,IWL,dn0,rtp,fthrd  &
  			    ,rtgt,f13t,f23t,glat,glon,rshort &
  			    ,rlong,albedt,cosz,rlongup,mynum  &
  			    ,fmapt,carma(ngrid)%aot,xland)
			    
!xxxxxxxxxxxxxx
!	    if (rlong(ii,jj).lt.10.0) print*,'apos radcarma',ii,jj,rlong(ii,jj)
!
!		 do k=2,m1-1
!		   if ( (fthrd(k,ii,jj))*86400 .lt. -20.0)then
!		    print*,'apos radcarma i j k mynum',ii,jj,k,mynum
!		    print*,'apos radcarma rlong fthrd',rlong(ii,jj),fthrd(k,ii,jj)*86400
!		   endif
!  	     end do
!xxxxxxxxxxxxxx
	    		    
  	 END DO
       END DO
       
       
       
  END SUBROUTINE radcomp_carma
  
    SUBROUTINE AllocIndex(ia,ja,iz,jz,IsAlloc)
      IMPLICIT NONE
      INTEGER,INTENT(IN) :: ia,iz,ja,jz,IsAlloc
      INTEGER :: ij,i1,j1
      
      IF(IsAlloc==1) THEN
  	ALLOCATE(indexi((iz-ia+1)*(jz-ja+1)))
  	ALLOCATE(indexj((iz-ia+1)*(jz-ja+1)))
      
  	ij=0
  	DO i1=ia,iz
  	  DO j1=ja,jz
  	    ij=ij+1
  	    indexi(ij)=i1
  	    indexj(ij)=j1
  	  END DO
  	END DO
      ELSE
  	DEALLOCATE(indexi)
  	DEALLOCATE(indexj)
      END IF
    
    END SUBROUTINE AllocIndex
    
    SUBROUTINE C_2d_1d(A2d,A1d,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend
       REAL,INTENT(IN) :: A2d(ia:iz,ja:jz)
       REAL,INTENT(OUT) :: A1d(iend)
       INTEGER :: ij
       
       DO ij=1,iend
  	 A1d(ij)=A2d(indexi(ij),indexj(ij))
       END DO
       
    END SUBROUTINE C_2d_1d
    
    SUBROUTINE C_1d_2d(A1d,A2d,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend
       REAL,INTENT(IN) :: A1d(iend)
       REAL,INTENT(OUT) :: A2d(ia:iz,ja:jz)
       INTEGER :: ij
       
       DO ij=1,iend
  	A2d(indexi(ij),indexj(ij))=A1d(ij)
       END DO
       
    END SUBROUTINE C_1d_2d
     
    SUBROUTINE C_3d_2d(A3d,A2d,m,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend,m
       REAL,INTENT(IN) :: A3d(m,ia:iz,ja:jz)
       REAL,INTENT(OUT) :: A2d(iend,m)
       INTEGER :: ij,l
       
       DO l=1,m
  	 DO ij=1,iend
  	   A2d(ij,l)=A3d(l,indexi(ij),indexj(ij))
  	 END DO
       END DO
       
    END SUBROUTINE C_3d_2d
  
    SUBROUTINE C_2d_3d(A2d,A3d,m,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend,m
       REAL,INTENT(OUT)   :: A3d(m,ia:iz,ja:jz)
       REAL,INTENT(IN)    :: A2d(iend,m)
       INTEGER :: ij,l
       
       DO l=1,m
  	 DO ij=1,iend
  	   A3d(l,indexi(ij),indexj(ij))=A2d(ij,l)
  	 END DO
       END DO
       
    END SUBROUTINE C_2d_3d
    
    SUBROUTINE C_3d_4d(A3d,A4d,m,ia,iz,ja,jz,iend,nk)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend,m,nk
       REAL,INTENT(OUT)   :: A4d(nk,m,ia:iz,ja:jz)
       REAL,INTENT(IN)    :: A3d(nk,iend,m)
       INTEGER :: ij,l,k
       
       DO l=1,m
  	 DO ij=1,iend
	   DO k=1,nk
  	     A4d(k,l,indexi(ij),indexj(ij))=A3d(k,ij,l)
	   END DO
	 END DO
       END DO
       
    END SUBROUTINE C_3d_4d
    
    
    SUBROUTINE Ci_3d_2d(A3d,A2d,m,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend,m
       REAL,INTENT(IN) :: A3d(ia:iz,ja:jz,m)
       REAL,INTENT(OUT) :: A2d(iend,m)
       INTEGER :: ij,l
       
       DO l=1,m
  	 DO ij=1,iend
  	   A2d(ij,l)=A3d(indexi(ij),indexj(ij),l)
  	 END DO
       END DO
       
    END SUBROUTINE Ci_3d_2d

    SUBROUTINE Ci_2d_3d(A2d,A3d,m,ia,iz,ja,jz,iend)
  
       INTEGER,INTENT(IN) :: ia,iz,ja,jz,iend,m
       REAL,INTENT(OUT)   :: A3d(ia:iz,ja:jz,m)
       REAL,INTENT(IN)    :: A2d(iend,m)
       INTEGER :: ij,l
       
       DO l=1,m
  	 DO ij=1,iend
  	   A3d(indexi(ij),indexj(ij),l)=A2d(ij,l)
  	 END DO
       END DO
       
    END SUBROUTINE Ci_2d_3d

  
END MODULE rad_carma
